Introduction

This markdown contains the full-length 16S rRNA amplicon sequencing analysis of 240 micro-ecosystems, which were treated with 4 stressors in all possible combinations at three different temperatures. The analysis is related to the manuscript “Stratified microbial communities are highly sensitive towards multiple combined global change factors, revealing antagonistic and synergistic effects” by Suleiman et al.

library(Biostrings)
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which, which.max, which.min
## Loading required package: S4Vectors
## Loading required package: stats4
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:base':
## 
##     expand.grid
## Loading required package: IRanges
## Loading required package: XVector
## 
## Attaching package: 'Biostrings'
## The following object is masked from 'package:base':
## 
##     strsplit
library(ShortRead)
## Loading required package: BiocParallel
## Loading required package: Rsamtools
## Loading required package: GenomeInfoDb
## Loading required package: GenomicRanges
## Loading required package: GenomicAlignments
## Loading required package: SummarizedExperiment
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: DelayedArray
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians
## 
## Attaching package: 'DelayedArray'
## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
## The following objects are masked from 'package:base':
## 
##     aperm, apply, rowsum
library(reshape2)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:Biobase':
## 
##     combine
## The following object is masked from 'package:BiocGenerics':
## 
##     combine
library(phyloseq)
## 
## Attaching package: 'phyloseq'
## The following object is masked from 'package:SummarizedExperiment':
## 
##     distance
## The following object is masked from 'package:Biobase':
## 
##     sampleNames
## The following object is masked from 'package:GenomicRanges':
## 
##     distance
## The following object is masked from 'package:IRanges':
## 
##     distance
library(boot)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.0     ✓ dplyr   1.0.4
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::collapse()   masks Biostrings::collapse(), IRanges::collapse()
## x dplyr::combine()    masks gridExtra::combine(), Biobase::combine(), BiocGenerics::combine()
## x purrr::compact()    masks XVector::compact()
## x purrr::compose()    masks ShortRead::compose()
## x dplyr::count()      masks matrixStats::count()
## x dplyr::desc()       masks IRanges::desc()
## x tidyr::expand()     masks S4Vectors::expand()
## x dplyr::filter()     masks stats::filter()
## x dplyr::first()      masks GenomicAlignments::first(), S4Vectors::first()
## x dplyr::id()         masks ShortRead::id()
## x dplyr::lag()        masks stats::lag()
## x dplyr::last()       masks GenomicAlignments::last()
## x ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## x purrr::reduce()     masks GenomicRanges::reduce(), IRanges::reduce()
## x dplyr::rename()     masks S4Vectors::rename()
## x purrr::simplify()   masks DelayedArray::simplify()
## x dplyr::slice()      masks XVector::slice(), IRanges::slice()
## x tibble::view()      masks ShortRead::view()
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:GenomicAlignments':
## 
##     second
## The following objects are masked from 'package:GenomicRanges':
## 
##     intersect, setdiff, union
## The following object is masked from 'package:GenomeInfoDb':
## 
##     intersect
## The following objects are masked from 'package:Biostrings':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:IRanges':
## 
##     %within%, intersect, setdiff, union
## The following objects are masked from 'package:S4Vectors':
## 
##     intersect, second, second<-, setdiff, union
## The following objects are masked from 'package:BiocGenerics':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(googlesheets)
library(here)
## here() starts at /Users/marcel/OneDrive - Universität Zürich UZH/Experiment FK/multiple_stressors/multiple_stressor
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
library(grid)
library(readxl)
library(ggplot2)
library(dplyr)
library(ggfortify)
library(ggpubr)
library(aod)
library(Rcpp)
library(dada2)
library(microbiome)
## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2020 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>
## 
## Attaching package: 'microbiome'
## The following object is masked from 'package:scales':
## 
##     alpha
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## The following object is masked from 'package:ShortRead':
## 
##     coverage
## The following object is masked from 'package:GenomicAlignments':
## 
##     coverage
## The following object is masked from 'package:SummarizedExperiment':
## 
##     coverage
## The following object is masked from 'package:GenomicRanges':
## 
##     coverage
## The following object is masked from 'package:Biostrings':
## 
##     coverage
## The following objects are masked from 'package:IRanges':
## 
##     coverage, transform
## The following object is masked from 'package:S4Vectors':
## 
##     transform
## The following object is masked from 'package:base':
## 
##     transform
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## 
## Attaching package: 'lattice'
## The following object is masked from 'package:boot':
## 
##     melanoma
## This is vegan 2.5-7
## 
## Attaching package: 'vegan'
## The following object is masked from 'package:microbiome':
## 
##     diversity
library(plyr)
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following object is masked from 'package:ggpubr':
## 
##     mutate
## The following object is masked from 'package:here':
## 
##     here
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:purrr':
## 
##     compact
## The following object is masked from 'package:ShortRead':
## 
##     id
## The following object is masked from 'package:matrixStats':
## 
##     count
## The following object is masked from 'package:XVector':
## 
##     compact
## The following object is masked from 'package:IRanges':
## 
##     desc
## The following object is masked from 'package:S4Vectors':
## 
##     rename
library(DESeq2)
library(apeglm)
library(ggpubr)
library(patchwork)
library(broom)
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following objects are masked from 'package:plyr':
## 
##     desc, mutate
## The following object is masked from 'package:IRanges':
## 
##     desc
## The following object is masked from 'package:stats':
## 
##     filter
library(plyr)
library(permute)
library(lattice)
library(tm)
## Loading required package: NLP
## 
## Attaching package: 'NLP'
## The following object is masked from 'package:microbiome':
## 
##     meta
## The following object is masked from 'package:ggplot2':
## 
##     annotate
## The following objects are masked from 'package:Biobase':
## 
##     annotation, content
## The following object is masked from 'package:BiocGenerics':
## 
##     annotation
## 
## Attaching package: 'tm'
## The following object is masked from 'package:ShortRead':
## 
##     reader
library(stringr)
library(viridis)           
## Loading required package: viridisLite
## 
## Attaching package: 'viridis'
## The following object is masked from 'package:scales':
## 
##     viridis_pal

Load phyloseq data, remove chloroplast sequences, remove samples that did not work

ps <- readRDS("phyloseq_multiple_stressor.rds")
taxa_names(ps) <- paste0("Seq", seq(ntaxa(ps)))

ps_new <-  subset_taxa(ps, Order != "Chloroplast")

ps_new_1 <- subset_samples(ps_new,Lable != "ori_1")
ps_new_2 <- subset_samples(ps_new_1,Lable != "ori_4")
ps_new_3 <- subset_samples(ps_new_2,Lable != "B_1_24")
ps_new_4 <- subset_samples(ps_new_3,Lable != "C_1_20")
ps_new_5 <- subset_samples(ps_new_4,Lable != "BCD_5_24")
ps_new_6 <- subset_samples(ps_new_5,Lable != "CD_1_24")
ps_new_7 <- subset_samples(ps_new_6,Lable != "BE_5_24")
ps_new_8 <- subset_samples(ps_new_7,Lable != "BCDE_5_24")
ps_new_9 <- subset_samples(ps_new_8,Lable != "CD_5_24")
ps_new_10 <- subset_samples(ps_new_9,Lable != "BE_5_20")
ps_new_11 <- subset_samples(ps_new_10,Lable != "CE_1_24")

Calculation of alpha diversity index

diversity_test<- estimate_richness(ps_new_11, split=TRUE, measures=NULL)
## Warning in estimate_richness(ps_new_11, split = TRUE, measures = NULL): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.
diversity_test$File <- rownames(diversity_test)
diversity_test$File <- gsub("\\.\\.", "--", diversity_test$File)
div_raw <- left_join(diversity_test,sample_data(ps_new_11)) 
## Joining, by = "File"
## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 object
div_raw <-  div_raw %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))%>%
  mutate(Temperature_factor=Temperature_factor)

#filter samples with low amount of reads
div_raw <- div_raw %>%
  filter(Observed >=100) %>%
  filter(Combination != "ori") %>%
  mutate(T2=Temperature^2,T=Temperature_factor) %>%
  select(-Temperature)%>%
  reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))

stderr <- function(x, na.rm=FALSE) {
  if (na.rm) x <- na.omit(x)
  sqrt(var(x)/length(x))
}

#calculating mean of shannon and observed reads
div_means <- div_raw %>%
  filter(T %in% c("20", "24", "28")) %>%
    dplyr::group_by(Treatment, Combinations, T, Stressor,Stressor_numeric) %>%
    dplyr::summarize(Shannon_mean = mean(Shannon),
                     Observed_mean = mean(Observed),
                     Shannon_SE = stderr(Shannon),
                     Observed_SE = stderr(Observed))%>%
  reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
## `summarise()` has grouped output by 'Treatment', 'Combinations', 'T', 'Stressor'. You can override using the `.groups` argument.
# plots

plot_diversity_all <- div_raw %>%
  filter(T!="19")%>%
ggplot() +
  geom_point(aes(x=Combinations, y=Shannon, shape=Stressor,col=Combinations), size = 3,position = position_dodge(width = 0.4)) +
  facet_wrap("T")+
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold")) +
  ylab ("Shannon index") +
  xlab("Stressor(s) applied") +
   scale_colour_viridis_d(direction = -1)+
   theme(axis.text.x = element_text(angle = 50, hjust = 1))


plot_diversity_means <- div_means %>%
ggplot() +
  geom_point(aes(x=T, y=Shannon_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=Shannon_mean,col=Combinations,ymin=Shannon_mean-Shannon_SE, ymax=Shannon_mean+Shannon_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=T, y=Shannon_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature [°C]") +
  theme(legend.position = "none") +
  ylab("Diversity_richness (mean)")

plot_diversity_means_numeric <- div_means %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=Shannon_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=Shannon_mean,col=T,ymin=Shannon_mean-Shannon_SE, ymax=Shannon_mean+Shannon_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold")) +
  xlab("Number of stressors") +
  ylab("Diversity_richness (mean)") +
  theme(legend.position = "none")

#final plot
(plot_diversity_means +  plot_diversity_means_numeric)

#shannon on combinations
anova_diversity <- lm(Shannon ~F*G*M*A*T, data = div_raw)
autoplot(anova_diversity)
## Warning: `arrange_()` was deprecated in dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help

hist(residuals(anova_diversity ))

anova(anova_diversity )
## Analysis of Variance Table
## 
## Response: Shannon
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## F           1  4.155  4.1551 16.0171 9.316e-05 ***
## G           1  0.109  0.1093  0.4214   0.51709    
## M           1  0.009  0.0090  0.0347   0.85238    
## A           1  6.820  6.8197 26.2889 7.853e-07 ***
## T           2  3.539  1.7694  6.8208   0.00141 ** 
## F:G         1  0.200  0.2003  0.7721   0.38080    
## F:M         1  0.055  0.0551  0.2125   0.64538    
## G:M         1  0.076  0.0761  0.2935   0.58869    
## F:A         1  4.638  4.6376 17.8770 3.820e-05 ***
## G:A         1  0.278  0.2778  1.0710   0.30216    
## M:A         1  0.041  0.0408  0.1572   0.69228    
## F:T         2  0.571  0.2853  1.0997   0.33531    
## G:T         2  0.137  0.0684  0.2636   0.76857    
## M:T         2  0.159  0.0797  0.3072   0.73594    
## A:T         2  0.041  0.0204  0.0785   0.92455    
## F:G:M       1  0.131  0.1313  0.5060   0.47784    
## F:G:A       1  1.239  1.2386  4.7744   0.03024 *  
## F:M:A       1  0.450  0.4500  1.7346   0.18958    
## G:M:A       1  0.361  0.3614  1.3930   0.23953    
## F:G:T       2  0.303  0.1516  0.5844   0.55852    
## F:M:T       2  0.009  0.0044  0.0170   0.98311    
## G:M:T       2  0.110  0.0548  0.2111   0.80990    
## F:A:T       2  1.498  0.7491  2.8877   0.05841 .  
## G:A:T       2  0.123  0.0616  0.2375   0.78884    
## M:A:T       2  0.336  0.1678  0.6467   0.52502    
## F:G:M:A     1  0.009  0.0087  0.0336   0.85481    
## F:G:M:T     2  0.470  0.2349  0.9055   0.40628    
## F:G:A:T     2  0.077  0.0387  0.1491   0.86158    
## F:M:A:T     2  0.069  0.0347  0.1336   0.87502    
## G:M:A:T     2  0.029  0.0143  0.0550   0.94646    
## F:G:M:A:T   2  0.470  0.2352  0.9067   0.40577    
## Residuals 172 44.619  0.2594                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_diversity)
## 
## Call:
## lm(formula = Shannon ~ F * G * M * A * T, data = div_raw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.93534 -0.20822  0.05585  0.27054  1.42868 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.521269   0.227778  24.240   <2e-16 ***
## F            0.074100   0.322127   0.230   0.8183    
## G            0.364889   0.341668   1.068   0.2870    
## M           -0.040310   0.341668  -0.118   0.9062    
## A            0.747909   0.341668   2.189   0.0299 *  
## T24          0.146260   0.322127   0.454   0.6504    
## T28         -0.180864   0.322127  -0.561   0.5752    
## F:G         -1.044270   0.469577  -2.224   0.0275 *  
## F:M         -0.109273   0.469577  -0.233   0.8163    
## G:M         -0.062624   0.496432  -0.126   0.8998    
## F:A         -0.526997   0.483191  -1.091   0.2769    
## G:A         -0.587326   0.483191  -1.216   0.2258    
## M:A         -0.166522   0.483191  -0.345   0.7308    
## F:T24       -0.186145   0.469577  -0.396   0.6923    
## F:T28        0.050948   0.455557   0.112   0.9111    
## G:T24       -0.478917   0.483191  -0.991   0.3230    
## G:T28       -0.393542   0.469577  -0.838   0.4032    
## M:T24       -0.444855   0.469577  -0.947   0.3448    
## M:T28       -0.263084   0.469577  -0.560   0.5760    
## A:T24       -0.103835   0.483191  -0.215   0.8301    
## A:T28        0.034482   0.469577   0.073   0.9415    
## F:G:M        0.656103   0.673778   0.974   0.3315    
## F:G:A        1.196845   0.683335   1.751   0.0816 .  
## F:M:A        0.006661   0.673778   0.010   0.9921    
## G:M:A        0.176000   0.683335   0.258   0.7971    
## F:G:T24      1.390935   0.673778   2.064   0.0405 *  
## F:G:T28      0.795634   0.654244   1.216   0.2256    
## F:M:T24      0.852003   0.673778   1.265   0.2078    
## F:M:T28      0.509086   0.654244   0.778   0.4376    
## G:M:T24      0.837978   0.708192   1.183   0.2383    
## G:M:T28      0.348712   0.673778   0.518   0.6054    
## F:A:T24     -0.001744   0.692761  -0.003   0.9980    
## F:A:T28     -0.493711   0.664083  -0.743   0.4582    
## G:A:T24      0.571687   0.708192   0.807   0.4206    
## G:A:T28      0.268085   0.664083   0.404   0.6869    
## M:A:T24      0.759147   0.683335   1.111   0.2681    
## M:A:T28      0.528232   0.664083   0.795   0.4275    
## F:G:M:A     -0.787214   0.952866  -0.826   0.4099    
## F:G:M:T24   -1.818117   0.988499  -1.839   0.0676 .  
## F:G:M:T28   -0.728523   0.932223  -0.781   0.4356    
## F:G:A:T24   -1.259151   0.984115  -1.279   0.2025    
## F:G:A:T28   -0.454171   0.939155  -0.484   0.6293    
## F:M:A:T24   -0.650923   0.966382  -0.674   0.5015    
## F:M:A:T28   -0.361584   0.932223  -0.388   0.6986    
## G:M:A:T24   -1.157196   1.001535  -1.155   0.2495    
## G:M:A:T28   -0.282198   0.939155  -0.300   0.7642    
## F:G:M:A:T24  1.807982   1.394852   1.296   0.1966    
## F:G:M:A:T28  0.425072   1.318363   0.322   0.7475    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5093 on 172 degrees of freedom
## Multiple R-squared:  0.3727, Adjusted R-squared:  0.2013 
## F-statistic: 2.174 on 47 and 172 DF,  p-value: 0.0001588

Calculation rel abundance, just sequences that appear > 0.01 %, making dataframes for each column position

ps_rel <- transform_sample_counts(ps_new_11, function(x) x / sum(x) )

relab_threshold <- 0.001

ps_relab <- filter_taxa(ps_rel, function(x) !(sum(x < relab_threshold) == length(x)), TRUE)
ntaxa(ps_new_11)
## [1] 20319
ntaxa(ps_relab)
## [1] 5830
ps_relative <- transform_sample_counts(ps_relab, function(x)  x / sum(x))

NMDS

mds_whole <- ps_relative@otu_table %>%
  as.data.frame() %>%
  metaMDS(., 
        distance = "bray",
        k = 3, ## number of dimensions to reduce to
        try = 100, ## number of random starts to try
        autotransform = FALSE ## best not to use
)
## Run 0 stress 0.1147651 
## Run 1 stress 0.117328 
## Run 2 stress 0.1201922 
## Run 3 stress 0.1180983 
## Run 4 stress 0.1171836 
## Run 5 stress 0.11507 
## ... Procrustes: rmse 0.006711343  max resid 0.08903122 
## Run 6 stress 0.1192138 
## Run 7 stress 0.1166306 
## Run 8 stress 0.1223326 
## Run 9 stress 0.1203552 
## Run 10 stress 0.1148027 
## ... Procrustes: rmse 0.002132423  max resid 0.01589763 
## Run 11 stress 0.1158032 
## Run 12 stress 0.1211516 
## Run 13 stress 0.1147548 
## ... New best solution
## ... Procrustes: rmse 0.002275816  max resid 0.02983052 
## Run 14 stress 0.1165917 
## Run 15 stress 0.1157718 
## Run 16 stress 0.1221672 
## Run 17 stress 0.1150052 
## ... Procrustes: rmse 0.005425105  max resid 0.06951436 
## Run 18 stress 0.1219735 
## Run 19 stress 0.1149344 
## ... Procrustes: rmse 0.00521367  max resid 0.0589545 
## Run 20 stress 0.1157971 
## Run 21 stress 0.1147687 
## ... Procrustes: rmse 0.002162335  max resid 0.02918938 
## Run 22 stress 0.1147639 
## ... Procrustes: rmse 0.002125425  max resid 0.02817898 
## Run 23 stress 0.1217448 
## Run 24 stress 0.1162657 
## Run 25 stress 0.1150722 
## ... Procrustes: rmse 0.007000091  max resid 0.08936665 
## Run 26 stress 0.1152504 
## ... Procrustes: rmse 0.006393371  max resid 0.07191682 
## Run 27 stress 0.1147718 
## ... Procrustes: rmse 0.002082273  max resid 0.0268729 
## Run 28 stress 0.1147541 
## ... New best solution
## ... Procrustes: rmse 0.0002902666  max resid 0.003265375 
## ... Similar to previous best
## Run 29 stress 0.1151349 
## ... Procrustes: rmse 0.006531455  max resid 0.06005487 
## Run 30 stress 0.1166066 
## Run 31 stress 0.1147726 
## ... Procrustes: rmse 0.001761228  max resid 0.02059767 
## Run 32 stress 0.1157936 
## Run 33 stress 0.1148025 
## ... Procrustes: rmse 0.004203744  max resid 0.05364829 
## Run 34 stress 0.1166647 
## Run 35 stress 0.1178293 
## Run 36 stress 0.1151829 
## ... Procrustes: rmse 0.006275269  max resid 0.06002028 
## Run 37 stress 0.1150307 
## ... Procrustes: rmse 0.00493667  max resid 0.04736456 
## Run 38 stress 0.1165166 
## Run 39 stress 0.1149577 
## ... Procrustes: rmse 0.00367769  max resid 0.03573953 
## Run 40 stress 0.1186977 
## Run 41 stress 0.1237333 
## Run 42 stress 0.1164721 
## Run 43 stress 0.1150399 
## ... Procrustes: rmse 0.006451346  max resid 0.08959954 
## Run 44 stress 0.1272646 
## Run 45 stress 0.1271736 
## Run 46 stress 0.1148437 
## ... Procrustes: rmse 0.003686582  max resid 0.04045134 
## Run 47 stress 0.1149457 
## ... Procrustes: rmse 0.004186228  max resid 0.04788836 
## Run 48 stress 0.1215899 
## Run 49 stress 0.1216396 
## Run 50 stress 0.1161779 
## Run 51 stress 0.1147746 
## ... Procrustes: rmse 0.002987398  max resid 0.04055853 
## Run 52 stress 0.1233454 
## Run 53 stress 0.1232819 
## Run 54 stress 0.1160627 
## Run 55 stress 0.1148072 
## ... Procrustes: rmse 0.002016763  max resid 0.02421261 
## Run 56 stress 0.1147672 
## ... Procrustes: rmse 0.001491363  max resid 0.01096601 
## Run 57 stress 0.1149019 
## ... Procrustes: rmse 0.004133279  max resid 0.03957484 
## Run 58 stress 0.1148842 
## ... Procrustes: rmse 0.003865038  max resid 0.03467656 
## Run 59 stress 0.1147631 
## ... Procrustes: rmse 0.00103653  max resid 0.007984618 
## ... Similar to previous best
## Run 60 stress 0.1158119 
## Run 61 stress 0.120041 
## Run 62 stress 0.1216382 
## Run 63 stress 0.1151596 
## ... Procrustes: rmse 0.006141564  max resid 0.08723947 
## Run 64 stress 0.1149327 
## ... Procrustes: rmse 0.00426985  max resid 0.05363242 
## Run 65 stress 0.12263 
## Run 66 stress 0.1227652 
## Run 67 stress 0.1219891 
## Run 68 stress 0.1163769 
## Run 69 stress 0.114763 
## ... Procrustes: rmse 0.00106261  max resid 0.0112492 
## Run 70 stress 0.1149026 
## ... Procrustes: rmse 0.005905025  max resid 0.05642366 
## Run 71 stress 0.1190239 
## Run 72 stress 0.116873 
## Run 73 stress 0.11577 
## Run 74 stress 0.1206862 
## Run 75 stress 0.1196646 
## Run 76 stress 0.1147574 
## ... Procrustes: rmse 0.0006005219  max resid 0.004559622 
## ... Similar to previous best
## Run 77 stress 0.1168291 
## Run 78 stress 0.1224913 
## Run 79 stress 0.1161734 
## Run 80 stress 0.1227812 
## Run 81 stress 0.1164785 
## Run 82 stress 0.1241439 
## Run 83 stress 0.1165355 
## Run 84 stress 0.1155676 
## Run 85 stress 0.1158396 
## Run 86 stress 0.1213225 
## Run 87 stress 0.1148268 
## ... Procrustes: rmse 0.002679179  max resid 0.02778069 
## Run 88 stress 0.1148036 
## ... Procrustes: rmse 0.002636979  max resid 0.03015428 
## Run 89 stress 0.1196198 
## Run 90 stress 0.1157682 
## Run 91 stress 0.1254822 
## Run 92 stress 0.1179103 
## Run 93 stress 0.1147765 
## ... Procrustes: rmse 0.001213848  max resid 0.00936875 
## ... Similar to previous best
## Run 94 stress 0.1148853 
## ... Procrustes: rmse 0.004558338  max resid 0.0541002 
## Run 95 stress 0.1169937 
## Run 96 stress 0.1230001 
## Run 97 stress 0.1201933 
## Run 98 stress 0.1168677 
## Run 99 stress 0.1221757 
## Run 100 stress 0.1271232 
## *** Solution reached
 ## 0.11

mds_whole_res <- ps_relative@sam_data %>%
  as.tibble() %>%
  select(Treatment, Stressor,Stressor_numeric, Temperature, Combinations, Temperature_factor, Lable, Name,F,G,M,A) %>%
  bind_cols(as.tibble(scores(mds_whole, display = "sites"))) %>%
  mutate(T2=Temperature^2,T=Temperature_factor) %>%
  select(-Temperature)
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 object
mds_whole_res <- mds_whole_res %>%
  filter(Temperature_factor != "19") %>%
  reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))

ggplot(mds_whole_res, aes(x = NMDS1, y = NMDS2)) +
  geom_point(size = 2) +
  facet_grid(Combinations~Temperature_factor)+
   theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))

Plots of NMDS distances that show high significance

#means of NMDS
means_NMDS <- mds_whole_res %>%
  filter(Temperature_factor %in% c("20", "24", "28")) %>%
    dplyr::group_by(Treatment,Combinations, Stressor, Stressor_numeric, Temperature_factor,F,G,M,A) %>%
    dplyr::summarize(NMDS1_mean = mean(NMDS1),
                     NMDS2_mean = mean(NMDS2),
                     NMDS3_mean = mean(NMDS3),
                      NMDS1_SE = stderr(NMDS1),
                     NMDS2_SE = stderr(NMDS2),
                     NMDS3_SE= stderr(NMDS3))
## `summarise()` has grouped output by 'Treatment', 'Combinations', 'Stressor', 'Stressor_numeric', 'Temperature_factor', 'F', 'G', 'M'. You can override using the `.groups` argument.

NMDS1 and NMDS2 plots

#NMDS1 supplement plot

nmds1_plot_temperature <- means_NMDS %>%
ggplot() +
  geom_point(aes(x=Temperature_factor, y=NMDS1_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Temperature_factor,y=NMDS1_mean,col=Combinations,ymin=NMDS1_mean-NMDS1_SE, ymax=NMDS1_mean+NMDS1_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=Temperature_factor, y=NMDS1_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature")  +
   ylab("NMDS1 (mean)") +
  
  theme(legend.position = "none")


plot_NMDS1_means_numeric <- means_NMDS %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=NMDS1_mean, col=Temperature_factor, shape= Temperature_factor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=NMDS1_mean,col=Temperature_factor,ymin=NMDS1_mean-NMDS1_SE, ymax=NMDS1_mean+NMDS1_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+  
  xlab("Number of stressors") +
  ylab("NMDS1 (mean)") +
  theme(legend.position = "none")


# plot NMDS1
nmds1_plot_temperature +  plot_NMDS1_means_numeric 

#NMDS2 supplement plot

NMDS2_plot_temperature <- means_NMDS %>%
ggplot() +
  geom_point(aes(x=Temperature_factor, y=NMDS2_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Temperature_factor,y=NMDS2_mean,col=Combinations,ymin=NMDS2_mean-NMDS2_SE, ymax=NMDS2_mean+NMDS2_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=Temperature_factor, y=NMDS2_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature") +
  theme(legend.position = "none") +
  ylab("NMDS2 (mean)")


plot_NMDS2_means_numeric <- means_NMDS %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=NMDS2_mean, col=Temperature_factor, shape= Temperature_factor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=NMDS2_mean,col=Temperature_factor,ymin=NMDS2_mean-NMDS2_SE, ymax=NMDS2_mean+NMDS2_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold")) +
  xlab("Number of stressors") +
  ylab("NMDS2 (mean)")+
  theme(legend.position = "none")

# plot NMDS2
NMDS2_plot_temperature +  plot_NMDS2_means_numeric 

###NMDS1 anova and plots

anova_NMDS1 <- lm(NMDS1 ~F*G*M*A*T, data = mds_whole_res)
autoplot(anova_NMDS1)

hist(residuals(anova_NMDS1 ))

anova(anova_NMDS1 )
## Analysis of Variance Table
## 
## Response: NMDS1
##            Df Sum Sq Mean Sq   F value    Pr(>F)    
## F           1 61.036  61.036 1028.3184 < 2.2e-16 ***
## G           1  0.002   0.002    0.0393  0.843076    
## M           1  0.029   0.029    0.4889  0.485297    
## A           1  2.625   2.625   44.2292 3.266e-10 ***
## T           2  0.652   0.326    5.4883  0.004844 ** 
## F:G         1  0.220   0.220    3.7031  0.055862 .  
## F:M         1  0.000   0.000    0.0009  0.976463    
## G:M         1  0.142   0.142    2.4002  0.123044    
## F:A         1  3.009   3.009   50.7008 2.375e-11 ***
## G:A         1  0.014   0.014    0.2341  0.629091    
## M:A         1  0.023   0.023    0.3956  0.530179    
## F:T         2  2.659   1.329   22.3970 1.991e-09 ***
## G:T         2  0.068   0.034    0.5768  0.562723    
## M:T         2  0.063   0.031    0.5270  0.591274    
## A:T         2  0.211   0.105    1.7773  0.172005    
## F:G:M       1  0.001   0.001    0.0203  0.886743    
## F:G:A       1  0.128   0.128    2.1539  0.143927    
## F:M:A       1  0.011   0.011    0.1884  0.664797    
## G:M:A       1  0.130   0.130    2.1874  0.140864    
## F:G:T       2  0.100   0.050    0.8402  0.433284    
## F:M:T       2  0.033   0.017    0.2821  0.754527    
## G:M:T       2  0.423   0.212    3.5659  0.030251 *  
## F:A:T       2  0.452   0.226    3.8116  0.023888 *  
## G:A:T       2  0.125   0.063    1.0551  0.350268    
## M:A:T       2  0.037   0.019    0.3119  0.732470    
## F:G:M:A     1  0.001   0.001    0.0243  0.876249    
## F:G:M:T     2  0.046   0.023    0.3859  0.680386    
## F:G:A:T     2  0.030   0.015    0.2526  0.777054    
## F:M:A:T     2  0.016   0.008    0.1310  0.877296    
## G:M:A:T     2  0.181   0.090    1.5243  0.220517    
## F:G:M:A:T   2  0.021   0.011    0.1779  0.837146    
## Residuals 183 10.862   0.059                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_NMDS1)
## 
## Call:
## lm(formula = NMDS1 ~ F * G * M * A * T, data = mds_whole_res)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90640 -0.12562 -0.00319  0.08877  0.98663 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.62357    0.10895  -5.723 4.21e-08 ***
## F            1.01436    0.15408   6.583 4.72e-10 ***
## G            0.08727    0.16343   0.534   0.5940    
## M           -0.10535    0.15408  -0.684   0.4950    
## A            0.25032    0.15408   1.625   0.1060    
## T24          0.04710    0.15408   0.306   0.7602    
## T28         -0.15322    0.15408  -0.994   0.3213    
## F:G         -0.03967    0.22462  -0.177   0.8600    
## F:M          0.15503    0.21791   0.711   0.4777    
## G:M         -0.07758    0.22462  -0.345   0.7302    
## F:A         -0.34377    0.22462  -1.531   0.1276    
## G:A         -0.08645    0.22462  -0.385   0.7008    
## M:A          0.07985    0.21791   0.366   0.7145    
## F:T24        0.11684    0.22462   0.520   0.6036    
## F:T28        0.51957    0.21791   2.384   0.0181 *  
## G:T24       -0.40315    0.22462  -1.795   0.0743 .  
## G:T28       -0.08559    0.22462  -0.381   0.7036    
## M:T24       -0.11534    0.21791  -0.529   0.5972    
## M:T28       -0.02001    0.21791  -0.092   0.9269    
## A:T24        0.30010    0.21791   1.377   0.1701    
## A:T28        0.20718    0.21791   0.951   0.3430    
## F:G:M       -0.03538    0.31295  -0.113   0.9101    
## F:G:A        0.32514    0.31765   1.024   0.3074    
## F:M:A       -0.12032    0.31295  -0.384   0.7011    
## G:M:A        0.01457    0.31295   0.047   0.9629    
## F:G:T24      0.19467    0.31765   0.613   0.5407    
## F:G:T28     -0.01761    0.31295  -0.056   0.9552    
## F:M:T24     -0.26511    0.31295  -0.847   0.3980    
## F:M:T28     -0.11238    0.30817  -0.365   0.7158    
## G:M:T24      0.62782    0.32535   1.930   0.0552 .  
## G:M:T28      0.16828    0.31295   0.538   0.5914    
## F:A:T24     -0.36116    0.32229  -1.121   0.2639    
## F:A:T28     -0.20663    0.31295  -0.660   0.5099    
## G:A:T24      0.18880    0.31765   0.594   0.5530    
## G:A:T28     -0.05402    0.31295  -0.173   0.8631    
## M:A:T24      0.20496    0.30817   0.665   0.5068    
## M:A:T28      0.02853    0.30817   0.093   0.9263    
## F:G:M:A     -0.07491    0.44257  -0.169   0.8658    
## F:G:M:T24    0.09541    0.45470   0.210   0.8340    
## F:G:M:T28    0.03307    0.43921   0.075   0.9401    
## F:G:A:T24   -0.33135    0.45252  -0.732   0.4650    
## F:G:A:T28   -0.01474    0.44257  -0.033   0.9735    
## F:M:A:T24    0.06044    0.44592   0.136   0.8923    
## F:M:A:T28    0.19036    0.43921   0.433   0.6652    
## G:M:A:T24   -0.62681    0.45143  -1.388   0.1667    
## G:M:A:T28    0.02288    0.43921   0.052   0.9585    
## F:G:M:A:T24  0.24337    0.63919   0.381   0.7038    
## F:G:M:A:T28 -0.13143    0.62114  -0.212   0.8327    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2436 on 183 degrees of freedom
## Multiple R-squared:  0.8697, Adjusted R-squared:  0.8362 
## F-statistic: 25.98 on 47 and 183 DF,  p-value: < 2.2e-16
#F:A:Temperature_factor  

ggplot(means_NMDS) + 
  geom_line(aes(x=Temperature_factor,y=NMDS1_mean, col=as.factor(F),shape=as.factor(A),
group=interaction(as.factor(F),A,M,G)), alpha=0.2)+
  geom_point(mapping = aes(x=Temperature_factor,y=NMDS1_mean, col=as.factor(F),shape=as.factor(A)),size=3)+
  theme_bw() + theme(panel.grid = element_blank())+
  geom_errorbar(aes(x=Temperature_factor,y=NMDS1_mean, ymin=NMDS1_mean-NMDS1_SE,ymax=NMDS1_mean+NMDS1_SE,col=as.factor(F)),width=0.1,alpha=0.2)+
  xlab("Temperature [°C]")
## Warning: Ignoring unknown aesthetics: shape

###NMDS2 anova and plots

#NMDS2 combination 
anova_NMDS2 <- lm(NMDS2 ~F*G*M*A*T, data = mds_whole_res)
autoplot(anova_NMDS2)

hist(residuals(anova_NMDS2 ))

anova(anova_NMDS2 )
## Analysis of Variance Table
## 
## Response: NMDS2
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## F           1  0.0493  0.0493  0.3960  0.529963    
## G           1  0.0004  0.0004  0.0034  0.953514    
## M           1  0.7552  0.7552  6.0673  0.014695 *  
## A           1  4.8989  4.8989 39.3559 2.484e-09 ***
## T           2  7.7550  3.8775 31.1500 2.275e-12 ***
## F:G         1  0.1971  0.1971  1.5831  0.209918    
## F:M         1  0.0523  0.0523  0.4200  0.517749    
## G:M         1  0.0059  0.0059  0.0471  0.828501    
## F:A         1  0.9378  0.9378  7.5338  0.006659 ** 
## G:A         1  0.4115  0.4115  3.3061  0.070658 .  
## M:A         1  0.1714  0.1714  1.3766  0.242201    
## F:T         2  0.0725  0.0362  0.2911  0.747764    
## G:T         2  0.2108  0.1054  0.8467  0.430497    
## M:T         2  0.2787  0.1394  1.1196  0.328627    
## A:T         2  0.2865  0.1433  1.1508  0.318652    
## F:G:M       1  0.0888  0.0888  0.7132  0.399501    
## F:G:A       1  0.0137  0.0137  0.1101  0.740426    
## F:M:A       1  0.3646  0.3646  2.9292  0.088686 .  
## G:M:A       1  0.4031  0.4031  3.2382  0.073587 .  
## F:G:T       2  0.0916  0.0458  0.3678  0.692757    
## F:M:T       2  0.1976  0.0988  0.7937  0.453731    
## G:M:T       2  0.3708  0.1854  1.4896  0.228179    
## F:A:T       2  1.1742  0.5871  4.7163  0.010064 *  
## G:A:T       2  0.4435  0.2217  1.7813  0.171333    
## M:A:T       2  0.1121  0.0561  0.4503  0.638123    
## F:G:M:A     1  0.0078  0.0078  0.0624  0.803066    
## F:G:M:T     2  0.0597  0.0299  0.2399  0.786919    
## F:G:A:T     2  0.2561  0.1281  1.0287  0.359528    
## F:M:A:T     2  0.0270  0.0135  0.1084  0.897316    
## G:M:A:T     2  0.1251  0.0625  0.5025  0.605855    
## F:G:M:A:T   2  0.7768  0.3884  3.1202  0.046506 *  
## Residuals 183 22.7795  0.1245                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_NMDS2)
## 
## Call:
## lm(formula = NMDS2 ~ F * G * M * A * T, data = mds_whole_res)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.67576 -0.10981 -0.00263  0.11760  1.33865 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.159311   0.157783  -1.010   0.3140  
## F            0.085881   0.223139   0.385   0.7008  
## G            0.067201   0.236675   0.284   0.7768  
## M            0.126427   0.223139   0.567   0.5717  
## A           -0.231951   0.223139  -1.039   0.2999  
## T24          0.434002   0.223139   1.945   0.0533 .
## T28          0.554193   0.223139   2.484   0.0139 *
## F:G          0.393576   0.325279   1.210   0.2279  
## F:M         -0.274407   0.315567  -0.870   0.3857  
## G:M         -0.153747   0.325279  -0.473   0.6370  
## F:A         -0.038816   0.325279  -0.119   0.9051  
## G:A          0.076456   0.325279   0.235   0.8144  
## M:A         -0.027077   0.315567  -0.086   0.9317  
## F:T24       -0.103381   0.325279  -0.318   0.7510  
## F:T28       -0.230081   0.315567  -0.729   0.4669  
## G:T24        0.144518   0.325279   0.444   0.6574  
## G:T28       -0.032296   0.325279  -0.099   0.9210  
## M:T24       -0.144553   0.315567  -0.458   0.6474  
## M:T28       -0.088370   0.315567  -0.280   0.7798  
## A:T24        0.094618   0.315567   0.300   0.7646  
## A:T28       -0.167913   0.315567  -0.532   0.5953  
## F:G:M       -0.335092   0.453198  -0.739   0.4606  
## F:G:A       -0.647511   0.460014  -1.408   0.1609  
## F:M:A        0.149664   0.453198   0.330   0.7416  
## G:M:A       -0.017552   0.453198  -0.039   0.9691  
## F:G:T24     -0.742368   0.460014  -1.614   0.1083  
## F:G:T28     -0.181871   0.453198  -0.401   0.6887  
## F:M:T24     -0.227877   0.453198  -0.503   0.6157  
## F:M:T28      0.192243   0.446279   0.431   0.6671  
## G:M:T24     -0.097784   0.471154  -0.208   0.8358  
## G:M:T28      0.031290   0.453198   0.069   0.9450  
## F:A:T24     -0.009746   0.466730  -0.021   0.9834  
## F:A:T28      0.541322   0.453198   1.194   0.2338  
## G:A:T24     -1.011734   0.460014  -2.199   0.0291 *
## G:A:T28     -0.090973   0.453198  -0.201   0.8411  
## M:A:T24     -0.449489   0.446279  -1.007   0.3152  
## M:A:T28     -0.044948   0.446279  -0.101   0.9199  
## F:G:M:A      0.580079   0.640919   0.905   0.3666  
## F:G:M:T24    1.183700   0.658482   1.798   0.0739 .
## F:G:M:T28    0.151148   0.636045   0.238   0.8124  
## F:G:A:T24    1.576498   0.655324   2.406   0.0171 *
## F:G:A:T28    0.138070   0.640919   0.215   0.8297  
## F:M:A:T24    0.691480   0.645756   1.071   0.2857  
## F:M:A:T28   -0.283331   0.636045  -0.445   0.6565  
## G:M:A:T24    0.766263   0.653739   1.172   0.2427  
## G:M:A:T28    0.153726   0.636045   0.242   0.8093  
## F:G:M:A:T24 -1.853952   0.925648  -2.003   0.0467 *
## F:G:M:A:T28  0.279283   0.899504   0.310   0.7565  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3528 on 183 degrees of freedom
## Multiple R-squared:  0.4748, Adjusted R-squared:  0.3399 
## F-statistic:  3.52 on 47 and 183 DF,  p-value: 6.71e-10

Visulaization of relative abundance of family level of all combinations with heatmap

df_rel <- psmelt(ps_relative)

df_family<- df_rel %>%
  group_by(Family,Lable,Combinations,Temperature_factor) %>%
  dplyr::summarize(Abundance=sum(Abundance))
## `summarise()` has grouped output by 'Family', 'Lable', 'Combinations'. You can override using the `.groups` argument.
df_family <- df_family %>%
  reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA")) 

df_family%>%
  filter(Family %in% c("Phormidiaceae","Chromatiaceae"),Combinations != "ori") %>%
  ggplot(aes(x=Combinations, y=Abundance, fill=Family)) +
  geom_boxplot()+
  facet_grid("Temperature_factor") +
  theme_bw()

df_family <- df_family %>%
  dplyr::group_by(Family,Combinations,Temperature_factor) %>%
  dplyr::summarize(Abundance_mean = mean(Abundance),
                   Abundance_SE =stderr(Abundance))
## `summarise()` has grouped output by 'Family', 'Combinations'. You can override using the `.groups` argument.
index <- which(df_family$Abundance_mean>=0.05)
family_to_keep <- unique(df_family[index,"Family"])
family_to_keep <- unname(unlist(family_to_keep))


df_family $Family_filter <- ifelse(df_family$Family %in% family_to_keep, df_family$Family,"other")

#make other to one
df_family<- df_family %>%
  dplyr::group_by(Family_filter,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance_sum = sum(Abundance_mean),
                  Abundance_sum_SE =sum(Abundance_SE))
## `summarise()` has grouped output by 'Family_filter', 'Combinations'. You can override using the `.groups` argument.
df_family <- df_family %>%
  reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA")) 

  cols <-c("Azospirillaceae"="bisque1", "Bacteroidetes_vadinHA17" = "yellow", "Chlorobiaceae"="black","Comamonadaceae" = "green4",  "Hungateiclostridiaceae" ="maroon1", "Lentimicrobiaceae" ="royalblue1","Marinifilaceae"="lightcyan1","Nocardiaceae" ="thistle1", "Phormidiaceae"= "firebrick1", "Prolixibacteraceae" = "magenta4", "Pseudomonadaceae" = "grey","Rhodocyclaceae"="green", "Rhodospillaceae"="grey38","Rubritaleaceae"="brown", "Xanthobacteraceae"="cornflowerblue", "Xanthomonadaceae"="lightskyblue1","Cyanobiaceae"="lavender", "Chromatiaceae" = "lightsalmon")

community_plot <- df_family %>%
  filter(Combinations != "ori",Family_filter !="Mitochondria",Family_filter !="other")%>%
ggplot(aes(x = Combinations, y = Abundance_sum, fill=Family_filter)) + 
    geom_col() +
  facet_wrap("Temperature_factor")+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 50, hjust = 1,size=14)) + 
    scale_fill_manual("Family", values = cols) +
  ylab("rel. abundance family level") +
  xlab("Stressor(s) applied")  +
  theme(axis.title=element_text(size=16,face="bold")) +
       theme(legend.text = element_text(size=11)) 

df_family %>%
  filter(Combinations != "ori", Family_filter != "other", Family_filter !="Pirellulaceae", Family_filter !="Cyanobiaceae",Family_filter !="Mitochondria")%>%
ggplot(aes(x = Combinations, y = Family_filter, fill=Abundance_sum)) + 
    geom_tile() +
    scale_fill_gradient2(low="red", mid="white", high="navy", midpoint=0.15, na.value="lightblue",
                 breaks=c(0.05,0.1,0.15,0.2,0.25),
                           limits=c(0.05,0.25)) +
  facet_wrap("Temperature_factor")+
  theme(strip.placement = "inside") +
  theme_bw()+
  theme(axis.text.x = element_text(angle = 50, hjust = 1))

#just phormidiaceae and chromatiaceae for supplementary information

df_family %>%
  filter(Family_filter %in% c("Phormidiaceae","Chromatiaceae"), Combinations!= "ori") %>%
  ggplot() +
  geom_point(aes(x=Combinations,y=Abundance_sum,col=Family_filter,fill=Family_filter),position=position_dodge(width = 0.1),size=4)+
  geom_errorbar(aes(x=Combinations,y=Abundance_sum,col=Family_filter,ymin=Abundance_sum-Abundance_sum_SE, ymax=Abundance_sum+Abundance_sum_SE),width=0.1,alpha=0.4,col="black", position=position_dodge(width = .1))+
  facet_wrap("Temperature_factor")+
  theme_bw() +
  ylab("rel. abundance") +
    theme(axis.text.x = element_text(angle = 50, hjust = 1, size=14)) + 
  labs(col = "Family") + labs(fill="Family")+
  xlab("Stressor combination")

Genus visualization

df_genus<- df_rel %>%
  group_by(Genus,Lable,Combinations,Temperature_factor) %>%
  dplyr::summarize(Abundance=sum(Abundance))
## `summarise()` has grouped output by 'Genus', 'Lable', 'Combinations'. You can override using the `.groups` argument.
df_genus <- df_genus %>%
  dplyr::group_by(Genus,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance = mean(Abundance))
## `summarise()` has grouped output by 'Genus', 'Combinations'. You can override using the `.groups` argument.
index <- which(df_genus$Abundance>=0.05)
genus_to_keep <- unique(df_genus[index,"Genus"])
genus_to_keep <- unname(unlist(genus_to_keep))

df_genus $Genus_filter <- ifelse(df_genus$Genus %in% genus_to_keep, df_genus$Genus,"other")

#make other to one
df_genus<- df_genus %>%
  dplyr::group_by(Genus_filter,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance = sum(Abundance))
## `summarise()` has grouped output by 'Genus_filter', 'Combinations'. You can override using the `.groups` argument.
df_genus <- df_genus %>%
  reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA")) 

genus_plot <- df_genus %>%
  filter(Combinations != "ori", Genus_filter != "other", Genus_filter !="Pirellulaceae", Genus_filter !="Cyanobium_PCC-6307", Genus_filter !="HN-HF0106")%>%
ggplot(aes(x = Combinations, y = Genus_filter, fill=Abundance)) + 
    geom_tile() +
  scale_fill_continuous(name = "rel.abundance",
                      low = "#FF0000",
                      high = "#00CCFF",
                 na.value="white",
                 breaks=c(0.05,0.1,0.15,0.2,0.25),
                           limits=c(0.05,0.25)) +
  facet_wrap("Temperature_factor") +
    theme(axis.text.x = element_text(angle = 50, hjust = 1, size=14)) +
   theme(axis.text.y = element_text(size=12,face="italic"))+ 
  theme(axis.title=element_text(size=16,face="bold")) +
      xlab("Stressor(s) applied") +
      ylab("rel. abundance genus level")
(community_plot) +
  (plot_spacer() + genus_plot + plot_layout(widths=c(1,100))) +
   plot_layout(ncol=1,heights=c(1,1)) + plot_annotation(tag_levels = "a")

Subpupulation dataset

ps_rel <- transform_sample_counts(ps_new_11, function(x) x / sum(x) )

relab_threshold <- 0.001

ps_relab <- filter_taxa(ps_rel, function(x) !(sum(x < relab_threshold) == length(x)), TRUE)
ntaxa(ps_new)
## [1] 20319
ntaxa(ps_relab)
## [1] 5830
ps_interest <- subset_taxa(ps_relab, Order %in% c("Cyanobacteriales","Chromatiales","Synechococcales", "Chlorobiales","Leptolyngbyales","Desulfobulbales","Desulfovibrionales","Limnotrichales","Campylobacterales"))

ntaxa(ps_interest)
## [1] 493
ps_relative_interest <- transform_sample_counts(ps_interest, function(x)  x / sum(x))

interest <-psmelt(ps_relative_interest)

ps_nmds <- subset_samples(ps_relative_interest, Treatment != "ori")

mds_whole_subset <- ps_nmds@otu_table %>%
  as.data.frame() %>%
  metaMDS(., 
        distance = "bray",
        k = 3, ## number of dimensions to reduce to
        try = 300, ## number of random starts to try
        autotransform = FALSE ## best not to use
) 
## Run 0 stress 0.1048781 
## Run 1 stress 0.1072036 
## Run 2 stress 0.1037414 
## ... New best solution
## ... Procrustes: rmse 0.006635281  max resid 0.08474109 
## Run 3 stress 0.105289 
## Run 4 stress 0.1051615 
## Run 5 stress 0.1061844 
## Run 6 stress 0.1063321 
## Run 7 stress 0.1072051 
## Run 8 stress 0.1049566 
## Run 9 stress 0.1052647 
## Run 10 stress 0.106756 
## Run 11 stress 0.1037896 
## ... Procrustes: rmse 0.00293663  max resid 0.02417046 
## Run 12 stress 0.1065792 
## Run 13 stress 0.1071902 
## Run 14 stress 0.1056136 
## Run 15 stress 0.1073591 
## Run 16 stress 0.1057231 
## Run 17 stress 0.1052825 
## Run 18 stress 0.1064098 
## Run 19 stress 0.1053836 
## Run 20 stress 0.1048949 
## Run 21 stress 0.1048907 
## Run 22 stress 0.1052388 
## Run 23 stress 0.1051941 
## Run 24 stress 0.1049397 
## Run 25 stress 0.105093 
## Run 26 stress 0.1037102 
## ... New best solution
## ... Procrustes: rmse 0.003714632  max resid 0.02826341 
## Run 27 stress 0.1054493 
## Run 28 stress 0.1053349 
## Run 29 stress 0.1037115 
## ... Procrustes: rmse 0.0005342914  max resid 0.007379712 
## ... Similar to previous best
## Run 30 stress 0.1039884 
## ... Procrustes: rmse 0.008394872  max resid 0.1217423 
## Run 31 stress 0.103687 
## ... New best solution
## ... Procrustes: rmse 0.001862688  max resid 0.02555876 
## Run 32 stress 0.1040644 
## ... Procrustes: rmse 0.008425929  max resid 0.1189921 
## Run 33 stress 0.1037651 
## ... Procrustes: rmse 0.002969157  max resid 0.04276497 
## Run 34 stress 0.1052573 
## Run 35 stress 0.1066791 
## Run 36 stress 0.1036805 
## ... New best solution
## ... Procrustes: rmse 0.0007859608  max resid 0.009224061 
## ... Similar to previous best
## Run 37 stress 0.1066239 
## Run 38 stress 0.1051974 
## Run 39 stress 0.1072624 
## Run 40 stress 0.103879 
## ... Procrustes: rmse 0.004549549  max resid 0.043179 
## Run 41 stress 0.1039829 
## ... Procrustes: rmse 0.007968812  max resid 0.117865 
## Run 42 stress 0.1054307 
## Run 43 stress 0.1065277 
## Run 44 stress 0.1063728 
## Run 45 stress 0.1071685 
## Run 46 stress 0.1053852 
## Run 47 stress 0.1051288 
## Run 48 stress 0.1072725 
## Run 49 stress 0.1055484 
## Run 50 stress 0.1054812 
## Run 51 stress 0.106358 
## Run 52 stress 0.1053356 
## Run 53 stress 0.1039857 
## ... Procrustes: rmse 0.008051974  max resid 0.1201322 
## Run 54 stress 0.1063092 
## Run 55 stress 0.1037752 
## ... Procrustes: rmse 0.002983071  max resid 0.04208786 
## Run 56 stress 0.1040888 
## ... Procrustes: rmse 0.008837648  max resid 0.1238711 
## Run 57 stress 0.1051795 
## Run 58 stress 0.1053434 
## Run 59 stress 0.1038061 
## ... Procrustes: rmse 0.003754823  max resid 0.0245486 
## Run 60 stress 0.1040709 
## ... Procrustes: rmse 0.007758431  max resid 0.1073494 
## Run 61 stress 0.1052072 
## Run 62 stress 0.107323 
## Run 63 stress 0.1051288 
## Run 64 stress 0.1037687 
## ... Procrustes: rmse 0.002974885  max resid 0.04227947 
## Run 65 stress 0.1037411 
## ... Procrustes: rmse 0.003008803  max resid 0.02138887 
## Run 66 stress 0.1067213 
## Run 67 stress 0.1052758 
## Run 68 stress 0.1052533 
## Run 69 stress 0.1096792 
## Run 70 stress 0.1054532 
## Run 71 stress 0.1051643 
## Run 72 stress 0.105183 
## Run 73 stress 0.1040851 
## ... Procrustes: rmse 0.008028575  max resid 0.1092318 
## Run 74 stress 0.1037936 
## ... Procrustes: rmse 0.003547532  max resid 0.04232409 
## Run 75 stress 0.1049463 
## Run 76 stress 0.1065428 
## Run 77 stress 0.1037843 
## ... Procrustes: rmse 0.003296965  max resid 0.04294894 
## Run 78 stress 0.1054323 
## Run 79 stress 0.1054745 
## Run 80 stress 0.1039897 
## ... Procrustes: rmse 0.007517483  max resid 0.1102222 
## Run 81 stress 0.1038034 
## ... Procrustes: rmse 0.003542826  max resid 0.02124059 
## Run 82 stress 0.1037411 
## ... Procrustes: rmse 0.00300768  max resid 0.02143947 
## Run 83 stress 0.1051659 
## Run 84 stress 0.1038814 
## ... Procrustes: rmse 0.004661334  max resid 0.04320149 
## Run 85 stress 0.1040868 
## ... Procrustes: rmse 0.008775891  max resid 0.1229213 
## Run 86 stress 0.1054308 
## Run 87 stress 0.1039924 
## ... Procrustes: rmse 0.008204912  max resid 0.122926 
## Run 88 stress 0.1036944 
## ... Procrustes: rmse 0.001268815  max resid 0.01526923 
## Run 89 stress 0.1052637 
## Run 90 stress 0.1055796 
## Run 91 stress 0.1052105 
## Run 92 stress 0.1040023 
## ... Procrustes: rmse 0.004446841  max resid 0.04310146 
## Run 93 stress 0.1040637 
## ... Procrustes: rmse 0.00846594  max resid 0.1184098 
## Run 94 stress 0.1039327 
## ... Procrustes: rmse 0.005352253  max resid 0.06325116 
## Run 95 stress 0.1044426 
## Run 96 stress 0.1052806 
## Run 97 stress 0.1054762 
## Run 98 stress 0.1065128 
## Run 99 stress 0.1052168 
## Run 100 stress 0.1086814 
## Run 101 stress 0.1052652 
## Run 102 stress 0.1053421 
## Run 103 stress 0.1055603 
## Run 104 stress 0.1053419 
## Run 105 stress 0.1049457 
## Run 106 stress 0.1066368 
## Run 107 stress 0.1064549 
## Run 108 stress 0.1037063 
## ... Procrustes: rmse 0.001999  max resid 0.02846922 
## Run 109 stress 0.1051921 
## Run 110 stress 0.1051284 
## Run 111 stress 0.103972 
## ... Procrustes: rmse 0.005570839  max resid 0.04323436 
## Run 112 stress 0.1054927 
## Run 113 stress 0.1054434 
## Run 114 stress 0.1056562 
## Run 115 stress 0.1066347 
## Run 116 stress 0.1057859 
## Run 117 stress 0.1051682 
## Run 118 stress 0.1066037 
## Run 119 stress 0.1053352 
## Run 120 stress 0.1054651 
## Run 121 stress 0.1055024 
## Run 122 stress 0.1040743 
## ... Procrustes: rmse 0.007981527  max resid 0.109647 
## Run 123 stress 0.1053739 
## Run 124 stress 0.1053455 
## Run 125 stress 0.106802 
## Run 126 stress 0.1039773 
## ... Procrustes: rmse 0.005982936  max resid 0.0510808 
## Run 127 stress 0.104097 
## ... Procrustes: rmse 0.008873997  max resid 0.1236601 
## Run 128 stress 0.1062863 
## Run 129 stress 0.104212 
## Run 130 stress 0.1053634 
## Run 131 stress 0.1071667 
## Run 132 stress 0.1037409 
## ... Procrustes: rmse 0.002918317  max resid 0.02096524 
## Run 133 stress 0.1051314 
## Run 134 stress 0.1068572 
## Run 135 stress 0.1051954 
## Run 136 stress 0.1053888 
## Run 137 stress 0.1069575 
## Run 138 stress 0.1054316 
## Run 139 stress 0.1082515 
## Run 140 stress 0.1053717 
## Run 141 stress 0.1065899 
## Run 142 stress 0.1039 
## ... Procrustes: rmse 0.005047809  max resid 0.06671632 
## Run 143 stress 0.1037197 
## ... Procrustes: rmse 0.002342273  max resid 0.02938035 
## Run 144 stress 0.103986 
## ... Procrustes: rmse 0.007920068  max resid 0.1183227 
## Run 145 stress 0.1051985 
## Run 146 stress 0.1039962 
## ... Procrustes: rmse 0.008331081  max resid 0.124856 
## Run 147 stress 0.1069593 
## Run 148 stress 0.1039855 
## ... Procrustes: rmse 0.008161781  max resid 0.1210844 
## Run 149 stress 0.1055484 
## Run 150 stress 0.1055342 
## Run 151 stress 0.103879 
## ... Procrustes: rmse 0.004546789  max resid 0.0431816 
## Run 152 stress 0.1087563 
## Run 153 stress 0.1072754 
## Run 154 stress 0.105523 
## Run 155 stress 0.1065537 
## Run 156 stress 0.1068461 
## Run 157 stress 0.1050111 
## Run 158 stress 0.1050574 
## Run 159 stress 0.1038874 
## ... Procrustes: rmse 0.004471144  max resid 0.04615632 
## Run 160 stress 0.1055009 
## Run 161 stress 0.1054307 
## Run 162 stress 0.1065876 
## Run 163 stress 0.1040893 
## ... Procrustes: rmse 0.008818907  max resid 0.1233156 
## Run 164 stress 0.1055367 
## Run 165 stress 0.1063708 
## Run 166 stress 0.1038614 
## ... Procrustes: rmse 0.005409715  max resid 0.07903686 
## Run 167 stress 0.104871 
## Run 168 stress 0.1055752 
## Run 169 stress 0.1038039 
## ... Procrustes: rmse 0.003589851  max resid 0.02210741 
## Run 170 stress 0.1039836 
## ... Procrustes: rmse 0.007939553  max resid 0.1172103 
## Run 171 stress 0.1037232 
## ... Procrustes: rmse 0.001996824  max resid 0.02116194 
## Run 172 stress 0.1053043 
## Run 173 stress 0.1051937 
## Run 174 stress 0.1049951 
## Run 175 stress 0.104861 
## Run 176 stress 0.1036858 
## ... Procrustes: rmse 0.0008393298  max resid 0.01193916 
## Run 177 stress 0.1064913 
## Run 178 stress 0.1038185 
## ... Procrustes: rmse 0.003540973  max resid 0.0421968 
## Run 179 stress 0.1044807 
## Run 180 stress 0.103804 
## ... Procrustes: rmse 0.003860474  max resid 0.0430969 
## Run 181 stress 0.1037412 
## ... Procrustes: rmse 0.003008108  max resid 0.02136756 
## Run 182 stress 0.1066865 
## Run 183 stress 0.1038053 
## ... Procrustes: rmse 0.003558565  max resid 0.04302693 
## Run 184 stress 0.1066373 
## Run 185 stress 0.1048932 
## Run 186 stress 0.1054485 
## Run 187 stress 0.1037931 
## ... Procrustes: rmse 0.003638536  max resid 0.04235056 
## Run 188 stress 0.105176 
## Run 189 stress 0.1049652 
## Run 190 stress 0.103821 
## ... Procrustes: rmse 0.004168878  max resid 0.04326228 
## Run 191 stress 0.1039914 
## ... Procrustes: rmse 0.006031939  max resid 0.07062848 
## Run 192 stress 0.1051735 
## Run 193 stress 0.1052762 
## Run 194 stress 0.1052536 
## Run 195 stress 0.1040649 
## ... Procrustes: rmse 0.008499819  max resid 0.1191021 
## Run 196 stress 0.1087645 
## Run 197 stress 0.106407 
## Run 198 stress 0.105535 
## Run 199 stress 0.1086678 
## Run 200 stress 0.1051626 
## Run 201 stress 0.1052206 
## Run 202 stress 0.1039534 
## ... Procrustes: rmse 0.006070404  max resid 0.07777913 
## Run 203 stress 0.1051645 
## Run 204 stress 0.1057991 
## Run 205 stress 0.1094136 
## Run 206 stress 0.1054887 
## Run 207 stress 0.1051517 
## Run 208 stress 0.1051306 
## Run 209 stress 0.1065712 
## Run 210 stress 0.1040688 
## ... Procrustes: rmse 0.008652522  max resid 0.1218283 
## Run 211 stress 0.1078543 
## Run 212 stress 0.1037077 
## ... Procrustes: rmse 0.002066305  max resid 0.02862835 
## Run 213 stress 0.1038549 
## ... Procrustes: rmse 0.004129181  max resid 0.04233274 
## Run 214 stress 0.1081731 
## Run 215 stress 0.1054627 
## Run 216 stress 0.1039943 
## ... Procrustes: rmse 0.008283253  max resid 0.1240307 
## Run 217 stress 0.1072255 
## Run 218 stress 0.1053806 
## Run 219 stress 0.1052059 
## Run 220 stress 0.1051759 
## Run 221 stress 0.1052599 
## Run 222 stress 0.1102512 
## Run 223 stress 0.1053471 
## Run 224 stress 0.1078586 
## Run 225 stress 0.105264 
## Run 226 stress 0.1051973 
## Run 227 stress 0.1042102 
## Run 228 stress 0.1063842 
## Run 229 stress 0.1040935 
## ... Procrustes: rmse 0.007944401  max resid 0.1085664 
## Run 230 stress 0.111074 
## Run 231 stress 0.1069745 
## Run 232 stress 0.1073528 
## Run 233 stress 0.103821 
## ... Procrustes: rmse 0.003117366  max resid 0.03421321 
## Run 234 stress 0.1055506 
## Run 235 stress 0.105613 
## Run 236 stress 0.1039835 
## ... Procrustes: rmse 0.007891517  max resid 0.1166033 
## Run 237 stress 0.1049836 
## Run 238 stress 0.1038503 
## ... Procrustes: rmse 0.003974036  max resid 0.04246804 
## Run 239 stress 0.1039039 
## ... Procrustes: rmse 0.005193599  max resid 0.04323721 
## Run 240 stress 0.1052643 
## Run 241 stress 0.1056454 
## Run 242 stress 0.1052456 
## Run 243 stress 0.1051273 
## Run 244 stress 0.103819 
## ... Procrustes: rmse 0.004125337  max resid 0.04318105 
## Run 245 stress 0.1039552 
## ... Procrustes: rmse 0.005854683  max resid 0.06957821 
## Run 246 stress 0.1054437 
## Run 247 stress 0.1058083 
## Run 248 stress 0.1053302 
## Run 249 stress 0.1052671 
## Run 250 stress 0.1070042 
## Run 251 stress 0.1039844 
## ... Procrustes: rmse 0.008057608  max resid 0.1196643 
## Run 252 stress 0.1037111 
## ... Procrustes: rmse 0.002050762  max resid 0.02656364 
## Run 253 stress 0.1055361 
## Run 254 stress 0.1040701 
## ... Procrustes: rmse 0.008383744  max resid 0.116903 
## Run 255 stress 0.105128 
## Run 256 stress 0.1051727 
## Run 257 stress 0.1058831 
## Run 258 stress 0.105211 
## Run 259 stress 0.1052793 
## Run 260 stress 0.1038193 
## ... Procrustes: rmse 0.004146237  max resid 0.0431907 
## Run 261 stress 0.105438 
## Run 262 stress 0.1050296 
## Run 263 stress 0.1042594 
## Run 264 stress 0.104895 
## Run 265 stress 0.1101512 
## Run 266 stress 0.1039903 
## ... Procrustes: rmse 0.007480987  max resid 0.1095657 
## Run 267 stress 0.103961 
## ... Procrustes: rmse 0.005597371  max resid 0.04327574 
## Run 268 stress 0.1049214 
## Run 269 stress 0.1063384 
## Run 270 stress 0.1039939 
## ... Procrustes: rmse 0.00801165  max resid 0.1181229 
## Run 271 stress 0.1060832 
## Run 272 stress 0.1067616 
## Run 273 stress 0.1051971 
## Run 274 stress 0.1065053 
## Run 275 stress 0.105126 
## Run 276 stress 0.1083203 
## Run 277 stress 0.105196 
## Run 278 stress 0.1040693 
## ... Procrustes: rmse 0.008306505  max resid 0.1179624 
## Run 279 stress 0.1065383 
## Run 280 stress 0.1055357 
## Run 281 stress 0.1082019 
## Run 282 stress 0.1055367 
## Run 283 stress 0.1062985 
## Run 284 stress 0.1038073 
## ... Procrustes: rmse 0.003760729  max resid 0.05310769 
## Run 285 stress 0.1037982 
## ... Procrustes: rmse 0.003700509  max resid 0.04244458 
## Run 286 stress 0.1040814 
## ... Procrustes: rmse 0.008594948  max resid 0.1197492 
## Run 287 stress 0.1052303 
## Run 288 stress 0.103711 
## ... Procrustes: rmse 0.00228947  max resid 0.03068486 
## Run 289 stress 0.1065655 
## Run 290 stress 0.1069532 
## Run 291 stress 0.1056077 
## Run 292 stress 0.1040067 
## ... Procrustes: rmse 0.006302214  max resid 0.04339678 
## Run 293 stress 0.1038044 
## ... Procrustes: rmse 0.003888339  max resid 0.04310388 
## Run 294 stress 0.106776 
## Run 295 stress 0.1092293 
## Run 296 stress 0.1039849 
## ... Procrustes: rmse 0.008079487  max resid 0.1200314 
## Run 297 stress 0.1053559 
## Run 298 stress 0.1070089 
## Run 299 stress 0.1041694 
## ... Procrustes: rmse 0.009284797  max resid 0.1249693 
## Run 300 stress 0.1054362 
## *** Solution reached
## 0.1

mds_whole_res_subset <- ps_nmds@sam_data %>%
  as.tibble() %>%
  select(Treatment,Temperature, Temperature_factor, Stressor,Stressor_numeric, Combinations, Temperature_factor, Lable, Name,F,G,M,A) %>%
  bind_cols(as.tibble(scores(mds_whole_subset, display = "sites"))) %>%
  mutate(T2=Temperature^2,T=Temperature_factor) %>%
  select(-Temperature)
## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 object
mds_whole_res_subset <- mds_whole_res_subset %>%
  reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))

# NMDS subset anova combination
#NMDS1 anova combination
anova_NMDS1_mixed <- lm(NMDS1 ~F*G*M*A*T, data = mds_whole_res_subset)
autoplot(anova_NMDS1)

hist(residuals(anova_NMDS1 ))

anova(anova_NMDS1 )
## Analysis of Variance Table
## 
## Response: NMDS1
##            Df Sum Sq Mean Sq   F value    Pr(>F)    
## F           1 61.036  61.036 1028.3184 < 2.2e-16 ***
## G           1  0.002   0.002    0.0393  0.843076    
## M           1  0.029   0.029    0.4889  0.485297    
## A           1  2.625   2.625   44.2292 3.266e-10 ***
## T           2  0.652   0.326    5.4883  0.004844 ** 
## F:G         1  0.220   0.220    3.7031  0.055862 .  
## F:M         1  0.000   0.000    0.0009  0.976463    
## G:M         1  0.142   0.142    2.4002  0.123044    
## F:A         1  3.009   3.009   50.7008 2.375e-11 ***
## G:A         1  0.014   0.014    0.2341  0.629091    
## M:A         1  0.023   0.023    0.3956  0.530179    
## F:T         2  2.659   1.329   22.3970 1.991e-09 ***
## G:T         2  0.068   0.034    0.5768  0.562723    
## M:T         2  0.063   0.031    0.5270  0.591274    
## A:T         2  0.211   0.105    1.7773  0.172005    
## F:G:M       1  0.001   0.001    0.0203  0.886743    
## F:G:A       1  0.128   0.128    2.1539  0.143927    
## F:M:A       1  0.011   0.011    0.1884  0.664797    
## G:M:A       1  0.130   0.130    2.1874  0.140864    
## F:G:T       2  0.100   0.050    0.8402  0.433284    
## F:M:T       2  0.033   0.017    0.2821  0.754527    
## G:M:T       2  0.423   0.212    3.5659  0.030251 *  
## F:A:T       2  0.452   0.226    3.8116  0.023888 *  
## G:A:T       2  0.125   0.063    1.0551  0.350268    
## M:A:T       2  0.037   0.019    0.3119  0.732470    
## F:G:M:A     1  0.001   0.001    0.0243  0.876249    
## F:G:M:T     2  0.046   0.023    0.3859  0.680386    
## F:G:A:T     2  0.030   0.015    0.2526  0.777054    
## F:M:A:T     2  0.016   0.008    0.1310  0.877296    
## G:M:A:T     2  0.181   0.090    1.5243  0.220517    
## F:G:M:A:T   2  0.021   0.011    0.1779  0.837146    
## Residuals 183 10.862   0.059                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_NMDS1)
## 
## Call:
## lm(formula = NMDS1 ~ F * G * M * A * T, data = mds_whole_res)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90640 -0.12562 -0.00319  0.08877  0.98663 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.62357    0.10895  -5.723 4.21e-08 ***
## F            1.01436    0.15408   6.583 4.72e-10 ***
## G            0.08727    0.16343   0.534   0.5940    
## M           -0.10535    0.15408  -0.684   0.4950    
## A            0.25032    0.15408   1.625   0.1060    
## T24          0.04710    0.15408   0.306   0.7602    
## T28         -0.15322    0.15408  -0.994   0.3213    
## F:G         -0.03967    0.22462  -0.177   0.8600    
## F:M          0.15503    0.21791   0.711   0.4777    
## G:M         -0.07758    0.22462  -0.345   0.7302    
## F:A         -0.34377    0.22462  -1.531   0.1276    
## G:A         -0.08645    0.22462  -0.385   0.7008    
## M:A          0.07985    0.21791   0.366   0.7145    
## F:T24        0.11684    0.22462   0.520   0.6036    
## F:T28        0.51957    0.21791   2.384   0.0181 *  
## G:T24       -0.40315    0.22462  -1.795   0.0743 .  
## G:T28       -0.08559    0.22462  -0.381   0.7036    
## M:T24       -0.11534    0.21791  -0.529   0.5972    
## M:T28       -0.02001    0.21791  -0.092   0.9269    
## A:T24        0.30010    0.21791   1.377   0.1701    
## A:T28        0.20718    0.21791   0.951   0.3430    
## F:G:M       -0.03538    0.31295  -0.113   0.9101    
## F:G:A        0.32514    0.31765   1.024   0.3074    
## F:M:A       -0.12032    0.31295  -0.384   0.7011    
## G:M:A        0.01457    0.31295   0.047   0.9629    
## F:G:T24      0.19467    0.31765   0.613   0.5407    
## F:G:T28     -0.01761    0.31295  -0.056   0.9552    
## F:M:T24     -0.26511    0.31295  -0.847   0.3980    
## F:M:T28     -0.11238    0.30817  -0.365   0.7158    
## G:M:T24      0.62782    0.32535   1.930   0.0552 .  
## G:M:T28      0.16828    0.31295   0.538   0.5914    
## F:A:T24     -0.36116    0.32229  -1.121   0.2639    
## F:A:T28     -0.20663    0.31295  -0.660   0.5099    
## G:A:T24      0.18880    0.31765   0.594   0.5530    
## G:A:T28     -0.05402    0.31295  -0.173   0.8631    
## M:A:T24      0.20496    0.30817   0.665   0.5068    
## M:A:T28      0.02853    0.30817   0.093   0.9263    
## F:G:M:A     -0.07491    0.44257  -0.169   0.8658    
## F:G:M:T24    0.09541    0.45470   0.210   0.8340    
## F:G:M:T28    0.03307    0.43921   0.075   0.9401    
## F:G:A:T24   -0.33135    0.45252  -0.732   0.4650    
## F:G:A:T28   -0.01474    0.44257  -0.033   0.9735    
## F:M:A:T24    0.06044    0.44592   0.136   0.8923    
## F:M:A:T28    0.19036    0.43921   0.433   0.6652    
## G:M:A:T24   -0.62681    0.45143  -1.388   0.1667    
## G:M:A:T28    0.02288    0.43921   0.052   0.9585    
## F:G:M:A:T24  0.24337    0.63919   0.381   0.7038    
## F:G:M:A:T28 -0.13143    0.62114  -0.212   0.8327    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2436 on 183 degrees of freedom
## Multiple R-squared:  0.8697, Adjusted R-squared:  0.8362 
## F-statistic: 25.98 on 47 and 183 DF,  p-value: < 2.2e-16
#NMDS2 anova combination
anova_NMDS2_mixed <- lm(NMDS2 ~F*G*M*A*T, data = mds_whole_res_subset)
autoplot(anova_NMDS2)

hist(residuals(anova_NMDS2 ))

anova(anova_NMDS2 )
## Analysis of Variance Table
## 
## Response: NMDS2
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## F           1  0.0493  0.0493  0.3960  0.529963    
## G           1  0.0004  0.0004  0.0034  0.953514    
## M           1  0.7552  0.7552  6.0673  0.014695 *  
## A           1  4.8989  4.8989 39.3559 2.484e-09 ***
## T           2  7.7550  3.8775 31.1500 2.275e-12 ***
## F:G         1  0.1971  0.1971  1.5831  0.209918    
## F:M         1  0.0523  0.0523  0.4200  0.517749    
## G:M         1  0.0059  0.0059  0.0471  0.828501    
## F:A         1  0.9378  0.9378  7.5338  0.006659 ** 
## G:A         1  0.4115  0.4115  3.3061  0.070658 .  
## M:A         1  0.1714  0.1714  1.3766  0.242201    
## F:T         2  0.0725  0.0362  0.2911  0.747764    
## G:T         2  0.2108  0.1054  0.8467  0.430497    
## M:T         2  0.2787  0.1394  1.1196  0.328627    
## A:T         2  0.2865  0.1433  1.1508  0.318652    
## F:G:M       1  0.0888  0.0888  0.7132  0.399501    
## F:G:A       1  0.0137  0.0137  0.1101  0.740426    
## F:M:A       1  0.3646  0.3646  2.9292  0.088686 .  
## G:M:A       1  0.4031  0.4031  3.2382  0.073587 .  
## F:G:T       2  0.0916  0.0458  0.3678  0.692757    
## F:M:T       2  0.1976  0.0988  0.7937  0.453731    
## G:M:T       2  0.3708  0.1854  1.4896  0.228179    
## F:A:T       2  1.1742  0.5871  4.7163  0.010064 *  
## G:A:T       2  0.4435  0.2217  1.7813  0.171333    
## M:A:T       2  0.1121  0.0561  0.4503  0.638123    
## F:G:M:A     1  0.0078  0.0078  0.0624  0.803066    
## F:G:M:T     2  0.0597  0.0299  0.2399  0.786919    
## F:G:A:T     2  0.2561  0.1281  1.0287  0.359528    
## F:M:A:T     2  0.0270  0.0135  0.1084  0.897316    
## G:M:A:T     2  0.1251  0.0625  0.5025  0.605855    
## F:G:M:A:T   2  0.7768  0.3884  3.1202  0.046506 *  
## Residuals 183 22.7795  0.1245                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_NMDS2)
## 
## Call:
## lm(formula = NMDS2 ~ F * G * M * A * T, data = mds_whole_res)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.67576 -0.10981 -0.00263  0.11760  1.33865 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.159311   0.157783  -1.010   0.3140  
## F            0.085881   0.223139   0.385   0.7008  
## G            0.067201   0.236675   0.284   0.7768  
## M            0.126427   0.223139   0.567   0.5717  
## A           -0.231951   0.223139  -1.039   0.2999  
## T24          0.434002   0.223139   1.945   0.0533 .
## T28          0.554193   0.223139   2.484   0.0139 *
## F:G          0.393576   0.325279   1.210   0.2279  
## F:M         -0.274407   0.315567  -0.870   0.3857  
## G:M         -0.153747   0.325279  -0.473   0.6370  
## F:A         -0.038816   0.325279  -0.119   0.9051  
## G:A          0.076456   0.325279   0.235   0.8144  
## M:A         -0.027077   0.315567  -0.086   0.9317  
## F:T24       -0.103381   0.325279  -0.318   0.7510  
## F:T28       -0.230081   0.315567  -0.729   0.4669  
## G:T24        0.144518   0.325279   0.444   0.6574  
## G:T28       -0.032296   0.325279  -0.099   0.9210  
## M:T24       -0.144553   0.315567  -0.458   0.6474  
## M:T28       -0.088370   0.315567  -0.280   0.7798  
## A:T24        0.094618   0.315567   0.300   0.7646  
## A:T28       -0.167913   0.315567  -0.532   0.5953  
## F:G:M       -0.335092   0.453198  -0.739   0.4606  
## F:G:A       -0.647511   0.460014  -1.408   0.1609  
## F:M:A        0.149664   0.453198   0.330   0.7416  
## G:M:A       -0.017552   0.453198  -0.039   0.9691  
## F:G:T24     -0.742368   0.460014  -1.614   0.1083  
## F:G:T28     -0.181871   0.453198  -0.401   0.6887  
## F:M:T24     -0.227877   0.453198  -0.503   0.6157  
## F:M:T28      0.192243   0.446279   0.431   0.6671  
## G:M:T24     -0.097784   0.471154  -0.208   0.8358  
## G:M:T28      0.031290   0.453198   0.069   0.9450  
## F:A:T24     -0.009746   0.466730  -0.021   0.9834  
## F:A:T28      0.541322   0.453198   1.194   0.2338  
## G:A:T24     -1.011734   0.460014  -2.199   0.0291 *
## G:A:T28     -0.090973   0.453198  -0.201   0.8411  
## M:A:T24     -0.449489   0.446279  -1.007   0.3152  
## M:A:T28     -0.044948   0.446279  -0.101   0.9199  
## F:G:M:A      0.580079   0.640919   0.905   0.3666  
## F:G:M:T24    1.183700   0.658482   1.798   0.0739 .
## F:G:M:T28    0.151148   0.636045   0.238   0.8124  
## F:G:A:T24    1.576498   0.655324   2.406   0.0171 *
## F:G:A:T28    0.138070   0.640919   0.215   0.8297  
## F:M:A:T24    0.691480   0.645756   1.071   0.2857  
## F:M:A:T28   -0.283331   0.636045  -0.445   0.6565  
## G:M:A:T24    0.766263   0.653739   1.172   0.2427  
## G:M:A:T28    0.153726   0.636045   0.242   0.8093  
## F:G:M:A:T24 -1.853952   0.925648  -2.003   0.0467 *
## F:G:M:A:T28  0.279283   0.899504   0.310   0.7565  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3528 on 183 degrees of freedom
## Multiple R-squared:  0.4748, Adjusted R-squared:  0.3399 
## F-statistic:  3.52 on 47 and 183 DF,  p-value: 6.71e-10
#NMDS3 anova combination
anova_NMDS3 <- lm(NMDS3 ~F*G*M*A*T, data = mds_whole_res_subset)
autoplot(anova_NMDS3)

hist(residuals(anova_NMDS3 ))

anova(anova_NMDS3 )
## Analysis of Variance Table
## 
## Response: NMDS3
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## F           1  0.1660 0.16603  1.1221   0.29087    
## G           1  0.0096 0.00956  0.0646   0.79969    
## M           1  0.0000 0.00005  0.0003   0.98601    
## A           1  0.2031 0.20307  1.3724   0.24292    
## T           2  3.9849 1.99245 13.4656 3.501e-06 ***
## F:G         1  0.0296 0.02957  0.1998   0.65537    
## F:M         1  0.4154 0.41544  2.8076   0.09552 .  
## G:M         1  0.4513 0.45134  3.0503   0.08240 .  
## F:A         1  0.2174 0.21738  1.4691   0.22705    
## G:A         1  0.9208 0.92075  6.2227   0.01350 *  
## M:A         1  0.3730 0.37299  2.5208   0.11408    
## F:T         2  5.6642 2.83212 19.1404 2.831e-08 ***
## G:T         2  0.0494 0.02469  0.1669   0.84645    
## M:T         2  0.4793 0.23963  1.6195   0.20082    
## A:T         2  0.3679 0.18394  1.2432   0.29090    
## F:G:M       1  0.4262 0.42623  2.8806   0.09135 .  
## F:G:A       1  0.0048 0.00481  0.0325   0.85716    
## F:M:A       1  0.0165 0.01645  0.1112   0.73917    
## G:M:A       1  0.2820 0.28196  1.9056   0.16914    
## F:G:T       2  0.3142 0.15709  1.0617   0.34799    
## F:M:T       2  0.0512 0.02560  0.1730   0.84125    
## G:M:T       2  0.0206 0.01028  0.0695   0.93293    
## F:A:T       2  0.0221 0.01107  0.0748   0.92793    
## G:A:T       2  0.1762 0.08810  0.5954   0.55242    
## M:A:T       2  0.4557 0.22784  1.5398   0.21719    
## F:G:M:A     1  0.1968 0.19683  1.3302   0.25027    
## F:G:M:T     2  0.2757 0.13783  0.9315   0.39583    
## F:G:A:T     2  0.0409 0.02044  0.1382   0.87105    
## F:M:A:T     2  0.0524 0.02619  0.1770   0.83792    
## G:M:A:T     2  1.2895 0.64477  4.3575   0.01417 *  
## F:G:M:A:T   2  0.4549 0.22747  1.5373   0.21772    
## Residuals 183 27.0778 0.14797                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_NMDS3)
## 
## Call:
## lm(formula = NMDS3 ~ F * G * M * A * T, data = mds_whole_res_subset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.08527 -0.18913 -0.00264  0.14258  1.25117 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.42168    0.17203   2.451  0.01517 * 
## F           -0.43074    0.24328  -1.771  0.07830 . 
## G           -0.05715    0.25804  -0.221  0.82498   
## M           -0.35646    0.24328  -1.465  0.14458   
## A           -0.07423    0.24328  -0.305  0.76063   
## T24         -0.50528    0.24328  -2.077  0.03921 * 
## T28         -0.77280    0.24328  -3.177  0.00175 **
## F:G          0.32615    0.35464   0.920  0.35896   
## F:M          0.15911    0.34405   0.462  0.64429   
## G:M          0.75950    0.35464   2.142  0.03355 * 
## F:A         -0.02962    0.35464  -0.084  0.93354   
## G:A         -0.06466    0.35464  -0.182  0.85553   
## M:A          0.61925    0.34405   1.800  0.07353 . 
## F:T24        0.61182    0.35464   1.725  0.08618 . 
## F:T28        0.48326    0.34405   1.405  0.16183   
## G:T24       -0.10935    0.35464  -0.308  0.75817   
## G:T28        0.09818    0.35464   0.277  0.78222   
## M:T24        0.21892    0.34405   0.636  0.52538   
## M:T28        0.38624    0.34405   1.123  0.26307   
## A:T24        0.40884    0.34405   1.188  0.23625   
## A:T28        0.31402    0.34405   0.913  0.36260   
## F:G:M       -0.99117    0.49411  -2.006  0.04633 * 
## F:G:A       -0.36071    0.50154  -0.719  0.47293   
## F:M:A       -0.21900    0.49411  -0.443  0.65813   
## G:M:A       -0.69423    0.49411  -1.405  0.16171   
## F:G:T24     -0.02092    0.50154  -0.042  0.96677   
## F:G:T28      0.09487    0.49411   0.192  0.84796   
## F:M:T24     -0.24707    0.49411  -0.500  0.61765   
## F:M:T28      0.03690    0.48657   0.076  0.93962   
## G:M:T24     -0.47948    0.51369  -0.933  0.35183   
## G:M:T28     -0.84546    0.49411  -1.711  0.08876 . 
## F:A:T24     -0.34859    0.50886  -0.685  0.49419   
## F:A:T28      0.34959    0.49411   0.708  0.48015   
## G:A:T24     -0.27441    0.50154  -0.547  0.58496   
## G:A:T28     -0.27332    0.49411  -0.553  0.58084   
## M:A:T24     -0.86023    0.48657  -1.768  0.07873 . 
## M:A:T28     -0.63003    0.48657  -1.295  0.19700   
## F:G:M:A      0.83301    0.69878   1.192  0.23476   
## F:G:M:T24    0.76706    0.71792   1.068  0.28673   
## F:G:M:T28    0.50833    0.69346   0.733  0.46448   
## F:G:A:T24    0.65627    0.71448   0.919  0.35955   
## F:G:A:T28   -0.35146    0.69878  -0.503  0.61559   
## F:M:A:T24    0.59664    0.70405   0.847  0.39785   
## F:M:A:T28   -0.41183    0.69346  -0.594  0.55333   
## G:M:A:T24    0.96117    0.71275   1.349  0.17916   
## G:M:A:T28    1.25184    0.69346   1.805  0.07269 . 
## F:G:M:A:T24 -1.41279    1.00921  -1.400  0.16323   
## F:G:M:A:T28  0.22341    0.98070   0.228  0.82005   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3847 on 183 degrees of freedom
## Multiple R-squared:  0.3914, Adjusted R-squared:  0.235 
## F-statistic: 2.504 on 47 and 183 DF,  p-value: 7.316e-06
#NMDS subset anova numeric

#NMDS3
anova_NMDS3_numeric <- lm(NMDS3 ~ Stressor_numeric*T, data = mds_whole_res_subset)
autoplot(anova_NMDS3_numeric)

hist(residuals(anova_NMDS3_numeric ))

anova(anova_NMDS3_numeric )
## Analysis of Variance Table
## 
## Response: NMDS3
##                     Df Sum Sq Mean Sq F value    Pr(>F)    
## Stressor_numeric     1  0.004 0.00405  0.0236  0.878077    
## T                    2  3.979 1.98941 11.5879 1.622e-05 ***
## Stressor_numeric:T   2  1.878 0.93918  5.4705  0.004788 ** 
## Residuals          225 38.628 0.17168                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_NMDS3_numeric)
## 
## Call:
## lm(formula = NMDS3 ~ Stressor_numeric * T, data = mds_whole_res_subset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17425 -0.22735 -0.04243  0.17112  1.46893 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           0.38990    0.10491   3.716 0.000255 ***
## Stressor_numeric     -0.10677    0.04662  -2.290 0.022936 *  
## T24                  -0.43546    0.14964  -2.910 0.003977 ** 
## T28                  -0.74761    0.14743  -5.071 8.27e-07 ***
## Stressor_numeric:T24  0.11091    0.06705   1.654 0.099482 .  
## Stressor_numeric:T28  0.21739    0.06572   3.308 0.001095 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4143 on 225 degrees of freedom
## Multiple R-squared:  0.1317, Adjusted R-squared:  0.1125 
## F-statistic: 6.828 on 5 and 225 DF,  p-value: 5.956e-06
#means of NMDS

means_NMDS_subset <- mds_whole_res_subset %>%
  filter(Temperature_factor %in% c("20", "24", "28")) %>%
    dplyr::group_by(Treatment,Combinations, Stressor, Stressor_numeric, Temperature_factor,F,G,M,A) %>%
    dplyr::summarize(NMDS1_mean = mean(NMDS1),
                     NMDS2_mean = mean(NMDS2),
                     NMDS3_mean = mean(NMDS3),
                      NMDS1_SE = stderr(NMDS1),
                     NMDS2_SE = stderr(NMDS2),
                     NMDS3_SE= stderr(NMDS3))
## `summarise()` has grouped output by 'Treatment', 'Combinations', 'Stressor', 'Stressor_numeric', 'Temperature_factor', 'F', 'G', 'M'. You can override using the `.groups` argument.
#NMDS1 subpopulation

nmds1_plot_temperature_subset <- means_NMDS_subset %>%
ggplot() +
  geom_point(aes(x=Temperature_factor, y=NMDS1_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Temperature_factor,y=NMDS1_mean,col=Combinations,ymin=NMDS1_mean-NMDS1_SE, ymax=NMDS1_mean+NMDS1_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=Temperature_factor, y=NMDS1_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature") +
  ylab("NMDS1_mean subpopulation")+
  theme(legend.position = "none") +
  ylab("NMDS1_subpopulation (mean)")


plot_NMDS1_means_numeric_subset<- means_NMDS_subset %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=NMDS1_mean, col=Temperature_factor, shape= Temperature_factor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=NMDS1_mean,col=Temperature_factor,ymin=NMDS1_mean-NMDS1_SE, ymax=NMDS1_mean+NMDS1_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
  ylab("NMDS1_subpopulation (mean)") +
  xlab( "Number of stressors")+
  theme(legend.position = "none")

# plot 
nmds1_plot_temperature_subset +  plot_NMDS1_means_numeric_subset 

#NMDS2 subpopulation

NMDS2_plot_temperature_subset <- means_NMDS_subset %>%
ggplot() +
  geom_point(aes(x=Temperature_factor, y=NMDS2_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Temperature_factor,y=NMDS2_mean,col=Combinations,ymin=NMDS2_mean-NMDS2_SE, ymax=NMDS2_mean+NMDS2_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=Temperature_factor, y=NMDS2_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature") +
 theme(legend.position = "none") +
  ylab("NMDS2_subpopulation (mean)")


plot_NMDS2_means_numeric_subset <- means_NMDS_subset %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=NMDS2_mean, col=Temperature_factor, shape= Temperature_factor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=NMDS2_mean,col=Temperature_factor,ymin=NMDS2_mean-NMDS2_SE, ymax=NMDS2_mean+NMDS2_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
  ylab("NMDS2_subpopulation (mean)") +
  xlab( "Number of stressors")+
 theme(legend.position = "none")

#plot 
NMDS2_plot_temperature_subset +  plot_NMDS2_means_numeric_subset 

Dataset of single genera, namely Tychonema_CCAP_1459-11B, Sulfuricurvum and Chromatium, calculating NMDS for all, respectively

ps_relative_cyano <- subset_taxa(ps_relative, Genus ==
                                "Tychonema_CCAP_1459-11B")

df_cyano <- psmelt(ps_relative_cyano)

df_cyano <- df_cyano %>%
  mutate(T2=Temperature^2,T=Temperature_factor) %>%
  select(-Temperature) 

df_cyano <- df_cyano %>%
  filter(Treatment != "ori")

df_cyano <- df_cyano %>%
dplyr::group_by(Treatment,Stressor_numeric,Lable,Combinations, T,F,G,M,A) %>%
dplyr::summarize(Abundance = sum(Abundance)) %>%
  reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
## `summarise()` has grouped output by 'Treatment', 'Stressor_numeric', 'Lable', 'Combinations', 'T', 'F', 'G', 'M'. You can override using the `.groups` argument.
#anova on combination
anova_tychonea <- lm(logit(Abundance) ~F*G*M*A*T, data = df_cyano)
autoplot(anova_tychonea)

hist(residuals(anova_tychonea))

anova(anova_tychonea)
## Analysis of Variance Table
## 
## Response: logit(Abundance)
##            Df Sum Sq Mean Sq  F value    Pr(>F)    
## F           1 93.372  93.372 191.0631 < 2.2e-16 ***
## G           1  0.037   0.037   0.0750  0.784453    
## M           1  0.304   0.304   0.6214  0.431543    
## A           1 79.227  79.227 162.1187 < 2.2e-16 ***
## T           2  4.268   2.134   4.3663  0.014047 *  
## F:G         1  0.062   0.062   0.1274  0.721557    
## F:M         1  0.014   0.014   0.0277  0.868010    
## G:M         1  2.988   2.988   6.1141  0.014323 *  
## F:A         1 77.297  77.297 158.1699 < 2.2e-16 ***
## G:A         1  0.391   0.391   0.8007  0.372060    
## M:A         1  2.395   2.395   4.9007  0.028081 *  
## F:T         2  6.685   3.342   6.8392  0.001366 ** 
## G:T         2  2.332   1.166   2.3864  0.094817 .  
## M:T         2  2.128   1.064   2.1775  0.116249    
## A:T         2  1.412   0.706   1.4450  0.238412    
## F:G:M       1  0.655   0.655   1.3403  0.248492    
## F:G:A       1  0.042   0.042   0.0859  0.769814    
## F:M:A       1  0.042   0.042   0.0851  0.770875    
## G:M:A       1  0.037   0.037   0.0767  0.782162    
## F:G:T       2  2.220   1.110   2.2709  0.106125    
## F:M:T       2  0.728   0.364   0.7443  0.476484    
## G:M:T       2  2.575   1.288   2.6346  0.074468 .  
## F:A:T       2  0.104   0.052   0.1067  0.898876    
## G:A:T       2  0.305   0.152   0.3120  0.732335    
## M:A:T       2  0.976   0.488   0.9985  0.370426    
## F:G:M:A     1  0.474   0.474   0.9696  0.326080    
## F:G:M:T     2  0.139   0.070   0.1423  0.867456    
## F:G:A:T     2  0.202   0.101   0.2068  0.813350    
## F:M:A:T     2  0.338   0.169   0.3454  0.708425    
## G:M:A:T     2  1.333   0.667   1.3638  0.258261    
## F:G:M:A:T   2  0.842   0.421   0.8612  0.424347    
## Residuals 183 89.432   0.489                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_tychonea)
## 
## Call:
## lm(formula = logit(Abundance) ~ F * G * M * A * T, data = df_cyano)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.48855 -0.31281 -0.00915  0.29578  2.41578 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.10233    0.31263  -6.725 2.18e-10 ***
## F           -2.08698    0.44213  -4.720 4.67e-06 ***
## G           -0.38890    0.46895  -0.829 0.408013    
## M            0.59739    0.44213   1.351 0.178310    
## A           -2.19776    0.44213  -4.971 1.53e-06 ***
## T24         -0.57596    0.44213  -1.303 0.194319    
## T28          0.13664    0.44213   0.309 0.757635    
## F:G          0.04950    0.64451   0.077 0.938860    
## F:M         -0.52908    0.62527  -0.846 0.398566    
## G:M         -0.76128    0.64451  -1.181 0.239068    
## F:A          2.16353    0.64451   3.357 0.000959 ***
## G:A         -0.12884    0.64451  -0.200 0.841776    
## M:A         -1.05617    0.62527  -1.689 0.092893 .  
## F:T24        0.11824    0.64451   0.183 0.854638    
## F:T28       -0.41784    0.62527  -0.668 0.504815    
## G:T24        1.65914    0.64451   2.574 0.010836 *  
## G:T28        0.60803    0.64451   0.943 0.346723    
## M:T24        0.84276    0.62527   1.348 0.179378    
## M:T28       -0.26699    0.62527  -0.427 0.669885    
## A:T24        0.36839    0.62527   0.589 0.556470    
## A:T28        0.01890    0.62527   0.030 0.975913    
## F:G:M        1.19959    0.89797   1.336 0.183243    
## F:G:A        0.57908    0.91147   0.635 0.526009    
## F:M:A        0.78792    0.89797   0.877 0.381396    
## G:M:A        1.35218    0.89797   1.506 0.133837    
## F:G:T24     -0.91244    0.91147  -1.001 0.318119    
## F:G:T28     -0.50304    0.89797  -0.560 0.576027    
## F:M:T24     -0.28087    0.89797  -0.313 0.754805    
## F:M:T28      0.38004    0.88426   0.430 0.667859    
## G:M:T24     -1.09066    0.93355  -1.168 0.244208    
## G:M:T28      0.75242    0.89797   0.838 0.403174    
## F:A:T24      0.01639    0.92478   0.018 0.985883    
## F:A:T28     -0.08825    0.89797  -0.098 0.921815    
## G:A:T24     -0.40545    0.91147  -0.445 0.656968    
## G:A:T28      0.72491    0.89797   0.807 0.420557    
## M:A:T24      0.51831    0.88426   0.586 0.558497    
## M:A:T28      0.53227    0.88426   0.602 0.547959    
## F:G:M:A     -1.94803    1.26992  -1.534 0.126762    
## F:G:M:T24   -0.19651    1.30472  -0.151 0.880446    
## F:G:M:T28   -1.05818    1.26026  -0.840 0.402200    
## F:G:A:T24   -0.11091    1.29846  -0.085 0.932025    
## F:G:A:T28   -0.85559    1.26992  -0.674 0.501332    
## F:M:A:T24   -0.50726    1.27951  -0.396 0.692235    
## F:M:A:T28   -0.45012    1.26026  -0.357 0.721381    
## G:M:A:T24   -0.32512    1.29532  -0.251 0.802097    
## G:M:A:T28   -2.25891    1.26026  -1.792 0.074719 .  
## F:G:M:A:T24  1.31007    1.83408   0.714 0.475958    
## F:G:M:A:T28  2.33497    1.78228   1.310 0.191804    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6991 on 183 degrees of freedom
## Multiple R-squared:  0.7605, Adjusted R-squared:  0.6989 
## F-statistic: 12.36 on 47 and 183 DF,  p-value: < 2.2e-16
# Sulfuricurvum
ps_sulfuri <- subset_taxa(ps_relative, Genus == "Sulfuricurvum")

df_sulfuri <- psmelt(ps_sulfuri)

df_sulfuri <- df_sulfuri %>%
  mutate(T2=Temperature^2,T=Temperature_factor) %>%
  select(-Temperature)

df_sulfuri <- df_sulfuri %>%
  filter(Treatment!="ori") %>%
dplyr::group_by(Treatment,Stressor_numeric,Combinations, Lable,T,F,G,M,A) %>%
dplyr::summarize(Abundance = sum(Abundance))
## `summarise()` has grouped output by 'Treatment', 'Stressor_numeric', 'Combinations', 'Lable', 'T', 'F', 'G', 'M'. You can override using the `.groups` argument.
#anova on combinations
anova_sulfuri <- lm(asin(sqrt(Abundance)) ~F*G*M*A*T, data = df_sulfuri)
autoplot(anova_sulfuri)

hist(residuals(anova_sulfuri))

anova(anova_sulfuri)
## Analysis of Variance Table
## 
## Response: asin(sqrt(Abundance))
##            Df  Sum Sq  Mean Sq F value    Pr(>F)    
## F           1 0.03216 0.032159 11.0700 0.0010604 ** 
## G           1 0.04274 0.042741 14.7128 0.0001722 ***
## M           1 0.00026 0.000256  0.0880 0.7670137    
## A           1 0.03900 0.038996 13.4236 0.0003253 ***
## T           2 0.03474 0.017369  5.9790 0.0030527 ** 
## F:G         1 0.00654 0.006539  2.2509 0.1352607    
## F:M         1 0.00647 0.006475  2.2287 0.1371887    
## G:M         1 0.00015 0.000153  0.0527 0.8187109    
## F:A         1 0.06559 0.065589 22.5777 4.069e-06 ***
## G:A         1 0.00159 0.001586  0.5460 0.4608872    
## M:A         1 0.00258 0.002580  0.8883 0.3471903    
## F:T         2 0.02910 0.014550  5.0084 0.0076261 ** 
## G:T         2 0.00094 0.000471  0.1621 0.8504703    
## M:T         2 0.01638 0.008189  2.8188 0.0622686 .  
## A:T         2 0.03860 0.019301  6.6439 0.0016388 ** 
## F:G:M       1 0.00183 0.001828  0.6293 0.4286264    
## F:G:A       1 0.00713 0.007128  2.4535 0.1189889    
## F:M:A       1 0.00007 0.000071  0.0243 0.8763079    
## G:M:A       1 0.00813 0.008128  2.7978 0.0961047 .  
## F:G:T       2 0.00806 0.004032  1.3878 0.2522327    
## F:M:T       2 0.02019 0.010097  3.4757 0.0329978 *  
## G:M:T       2 0.00117 0.000583  0.2007 0.8183574    
## F:A:T       2 0.00100 0.000499  0.1718 0.8423117    
## G:A:T       2 0.00371 0.001857  0.6392 0.5289039    
## M:A:T       2 0.01492 0.007458  2.5673 0.0795012 .  
## F:G:M:A     1 0.00034 0.000339  0.1166 0.7331165    
## F:G:M:T     2 0.01309 0.006547  2.2538 0.1079070    
## F:G:A:T     2 0.00449 0.002247  0.7735 0.4629033    
## F:M:A:T     2 0.00230 0.001151  0.3961 0.6735055    
## G:M:A:T     2 0.00395 0.001973  0.6791 0.5083566    
## F:G:M:A:T   2 0.00374 0.001871  0.6440 0.5263906    
## Residuals 183 0.53162 0.002905                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(anova_sulfuri)
## 
## Call:
## lm(formula = asin(sqrt(Abundance)) ~ F * G * M * A * T, data = df_sulfuri)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.144494 -0.020506 -0.005465  0.015081  0.259963 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.0708255  0.0241041   2.938  0.00372 **
## F           -0.0067231  0.0340884  -0.197  0.84387   
## G            0.0074519  0.0361562   0.206  0.83694   
## M           -0.0294875  0.0340884  -0.865  0.38815   
## A           -0.0503197  0.0340884  -1.476  0.14162   
## T24         -0.0543332  0.0340884  -1.594  0.11269   
## T28         -0.0654043  0.0340884  -1.919  0.05658 . 
## F:G          0.0654081  0.0496919   1.316  0.18973   
## F:M          0.0608374  0.0482083   1.262  0.20857   
## G:M          0.0000473  0.0496919   0.001  0.99924   
## F:A         -0.0055650  0.0496919  -0.112  0.91095   
## G:A          0.0255688  0.0496919   0.515  0.60749   
## M:A          0.0384315  0.0482083   0.797  0.42637   
## F:T24        0.1123383  0.0496919   2.261  0.02496 * 
## F:T28        0.0187107  0.0482083   0.388  0.69838   
## G:T24        0.0041848  0.0496919   0.084  0.93298   
## G:T28        0.0224243  0.0496919   0.451  0.65233   
## M:T24        0.0524766  0.0482083   1.089  0.27779   
## M:T28        0.0495668  0.0482083   1.028  0.30522   
## A:T24        0.0956222  0.0482083   1.984  0.04880 * 
## A:T28        0.0620968  0.0482083   1.288  0.19934   
## F:G:M       -0.0321723  0.0692338  -0.465  0.64271   
## F:G:A       -0.1011814  0.0702750  -1.440  0.15163   
## F:M:A       -0.0779992  0.0692338  -1.127  0.26138   
## G:M:A       -0.0365612  0.0692338  -0.528  0.59808   
## F:G:T24     -0.0322423  0.0702750  -0.459  0.64692   
## F:G:T28      0.0118470  0.0692338   0.171  0.86432   
## F:M:T24     -0.1272817  0.0692338  -1.838  0.06762 . 
## F:M:T28     -0.0604214  0.0681767  -0.886  0.37665   
## G:M:T24     -0.0178443  0.0719768  -0.248  0.80448   
## G:M:T28     -0.0195382  0.0692338  -0.282  0.77810   
## F:A:T24     -0.0916130  0.0713010  -1.285  0.20046   
## F:A:T28     -0.0091811  0.0692338  -0.133  0.89465   
## G:A:T24     -0.0301866  0.0702750  -0.430  0.66803   
## G:A:T28     -0.0631096  0.0692338  -0.912  0.36321   
## M:A:T24     -0.0900777  0.0681767  -1.321  0.18807   
## M:A:T28     -0.0278244  0.0681767  -0.408  0.68366   
## F:G:M:A      0.1082929  0.0979114   1.106  0.27017   
## F:G:M:T24    0.0555469  0.1005944   0.552  0.58149   
## F:G:M:T28   -0.0500372  0.0971668  -0.515  0.60720   
## F:G:A:T24    0.1073489  0.1001119   1.072  0.28500   
## F:G:A:T28    0.0362914  0.0979114   0.371  0.71132   
## F:M:A:T24    0.1125290  0.0986504   1.141  0.25549   
## F:M:A:T28    0.0792571  0.0971668   0.816  0.41574   
## G:M:A:T24    0.0671390  0.0998698   0.672  0.50226   
## G:M:A:T28    0.1496609  0.0971668   1.540  0.12523   
## F:G:M:A:T24 -0.1123707  0.1414086  -0.795  0.42784   
## F:G:M:A:T28 -0.1507164  0.1374146  -1.097  0.27417   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0539 on 183 degrees of freedom
## Multiple R-squared:  0.436,  Adjusted R-squared:  0.2911 
## F-statistic:  3.01 on 47 and 183 DF,  p-value: 7.15e-08

Abiotic factors and oxygen

abiotic<- ps_relative@sam_data %>%
  as.tibble() %>%
  select(Treatment,Temperature, Combinations, Temperature_factor, Lable, Name,F,G,M,A, Oxygen, TN,TC, pH,Stressor_numeric, Stressor) %>%
  mutate(T2=Temperature^2,T=Temperature_factor) %>%
  select(-Temperature)
## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 object
pH_anova <- lm(pH ~F*G*M*A*T, data = abiotic)
autoplot(pH_anova)

hist(residuals(pH_anova))

anova(pH_anova)
## Analysis of Variance Table
## 
## Response: pH
##            Df Sum Sq Mean Sq   F value    Pr(>F)    
## F           1 357.64  357.64 3216.0251 < 2.2e-16 ***
## G           1   0.14    0.14    1.2477 0.2654646    
## M           1   0.10    0.10    0.8570 0.3557852    
## A           1   1.77    1.77   15.9178 9.559e-05 ***
## T           2   0.83    0.42    3.7372 0.0256572 *  
## F:G         1   0.03    0.03    0.2885 0.5918137    
## F:M         1   0.27    0.27    2.3868 0.1240949    
## G:M         1   0.07    0.07    0.6723 0.4133055    
## F:A         1   1.62    1.62   14.6075 0.0001813 ***
## G:A         1   0.15    0.15    1.3826 0.2411803    
## M:A         1   0.43    0.43    3.9011 0.0497571 *  
## F:T         2   2.97    1.48   13.3530 3.863e-06 ***
## G:T         2   0.10    0.05    0.4283 0.6522488    
## M:T         2   0.02    0.01    0.0783 0.9247215    
## A:T         2   0.06    0.03    0.2584 0.7725663    
## F:G:M       1   0.05    0.05    0.4493 0.5035135    
## F:G:A       1   0.03    0.03    0.3072 0.5800555    
## F:M:A       1   0.36    0.36    3.2776 0.0718725 .  
## G:M:A       1   0.03    0.03    0.2406 0.6243562    
## F:G:T       2   0.82    0.41    3.6997 0.0266006 *  
## F:M:T       2   0.05    0.02    0.2117 0.8093793    
## G:M:T       2   0.05    0.03    0.2255 0.7983206    
## F:A:T       2   0.09    0.04    0.3823 0.6828679    
## G:A:T       2   0.54    0.27    2.4332 0.0905890 .  
## M:A:T       2   0.63    0.32    2.8499 0.0604239 .  
## F:G:M:A     1   0.01    0.01    0.0937 0.7598506    
## F:G:M:T     2   0.09    0.04    0.3829 0.6824486    
## F:G:A:T     2   0.05    0.02    0.2229 0.8004075    
## F:M:A:T     2   0.57    0.29    2.5786 0.0786348 .  
## G:M:A:T     2   0.12    0.06    0.5275 0.5909969    
## F:G:M:A:T   2   0.10    0.05    0.4664 0.6280159    
## Residuals 183  20.35    0.11                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(pH_anova)
## 
## Call:
## lm(formula = pH ~ F * G * M * A * T, data = abiotic)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.434 -0.102 -0.002  0.125  2.126 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.5600     0.1491  57.398   <2e-16 ***
## F            -2.6420     0.2109 -12.527   <2e-16 ***
## G            -0.0375     0.2237  -0.168   0.8671    
## M             0.2660     0.2109   1.261   0.2088    
## A             0.1860     0.2109   0.882   0.3790    
## T24           0.0440     0.2109   0.209   0.8350    
## T28           0.2480     0.2109   1.176   0.2412    
## F:G           0.2315     0.3075   0.753   0.4524    
## F:M          -0.2920     0.2983  -0.979   0.3289    
## G:M           0.0555     0.3075   0.181   0.8569    
## F:A          -0.1465     0.3075  -0.477   0.6343    
## G:A          -0.2065     0.3075  -0.672   0.5026    
## M:A          -0.7380     0.2983  -2.474   0.0143 *  
## F:T24         0.1405     0.3075   0.457   0.6482    
## F:T28        -0.0900     0.2983  -0.302   0.7632    
## G:T24         0.2075     0.3075   0.675   0.5006    
## G:T28         0.2375     0.3075   0.772   0.4408    
## M:T24        -0.3900     0.2983  -1.308   0.1927    
## M:T28        -0.1560     0.2983  -0.523   0.6016    
## A:T24        -0.7500     0.2983  -2.515   0.0128 *  
## A:T28        -0.3512     0.2983  -1.177   0.2405    
## F:G:M         0.2425     0.4284   0.566   0.5720    
## F:G:A         0.0330     0.4348   0.076   0.9396    
## F:M:A         0.8865     0.4284   2.070   0.0399 *  
## G:M:A        -0.1315     0.4284  -0.307   0.7592    
## F:G:T24      -0.3040     0.4348  -0.699   0.4853    
## F:G:T28      -0.4455     0.4284  -1.040   0.2997    
## F:M:T24       0.6215     0.4284   1.451   0.1485    
## F:M:T28       0.0760     0.4218   0.180   0.8572    
## G:M:T24      -0.2255     0.4453  -0.506   0.6132    
## G:M:T28      -0.1135     0.4284  -0.265   0.7913    
## F:A:T24       0.8305     0.4411   1.883   0.0613 .  
## F:A:T28       0.3357     0.4284   0.784   0.4342    
## G:A:T24       0.4165     0.4348   0.958   0.3394    
## G:A:T28       0.2657     0.4284   0.620   0.5358    
## M:A:T24       0.8760     0.4218   2.077   0.0392 *  
## M:A:T28       0.5512     0.4218   1.307   0.1929    
## F:G:M:A      -0.3350     0.6058  -0.553   0.5809    
## F:G:M:T24    -0.3505     0.6224  -0.563   0.5740    
## F:G:M:T28    -0.1805     0.6012  -0.300   0.7643    
## F:G:A:T24    -0.3390     0.6194  -0.547   0.5848    
## F:G:A:T28    -0.2022     0.6058  -0.334   0.7389    
## F:M:A:T24    -1.1965     0.6104  -1.960   0.0515 .  
## F:M:A:T28    -0.6417     0.6012  -1.067   0.2872    
## G:M:A:T24     0.2235     0.6179   0.362   0.7180    
## G:M:A:T28    -0.1797     0.6012  -0.299   0.7653    
## F:G:M:A:T24   0.4510     0.8749   0.515   0.6068    
## F:G:M:A:T28   0.8202     0.8502   0.965   0.3360    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3335 on 183 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.9478, Adjusted R-squared:  0.9344 
## F-statistic: 70.75 on 47 and 183 DF,  p-value: < 2.2e-16
oxygen_anova <- lm(Oxygen ~F*G*M*A*T, data = abiotic)
autoplot(oxygen_anova)

hist(residuals(oxygen_anova))

anova(oxygen_anova)
## Analysis of Variance Table
## 
## Response: Oxygen
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## F           1  24.59  24.594  6.9403 0.0091497 ** 
## G           1  21.31  21.315  6.0148 0.0151244 *  
## M           1  47.88  47.878 13.5108 0.0003115 ***
## A           1 161.67 161.674 45.6227 1.845e-10 ***
## T           2 183.84  91.919 25.9385 1.209e-10 ***
## F:G         1  12.87  12.866  3.6306 0.0582953 .  
## F:M         1   6.04   6.037  1.7036 0.1934602    
## G:M         1   1.37   1.366  0.3856 0.5353968    
## F:A         1  52.70  52.697 14.8705 0.0001593 ***
## G:A         1   4.33   4.334  1.2231 0.2702083    
## M:A         1   0.23   0.231  0.0651 0.7989686    
## F:T         2  34.02  17.012  4.8006 0.0092891 ** 
## G:T         2   9.56   4.779  1.3485 0.2622094    
## M:T         2   3.51   1.756  0.4955 0.6100538    
## A:T         2  50.91  25.454  7.1828 0.0009930 ***
## F:G:M       1  20.32  20.323  5.7349 0.0176401 *  
## F:G:A       1   8.83   8.828  2.4913 0.1162046    
## F:M:A       1   0.10   0.099  0.0278 0.8676644    
## G:M:A       1   0.00   0.001  0.0003 0.9870868    
## F:G:T       2  20.75  10.377  2.9283 0.0559984 .  
## F:M:T       2  17.88   8.939  2.5225 0.0830454 .  
## G:M:T       2   2.61   1.304  0.3681 0.6925778    
## F:A:T       2   0.31   0.155  0.0437 0.9572609    
## G:A:T       2   1.60   0.802  0.2264 0.7976449    
## M:A:T       2  18.25   9.123  2.5745 0.0789470 .  
## F:G:M:A     1  15.10  15.105  4.2624 0.0403759 *  
## F:G:M:T     2   5.41   2.705  0.7634 0.4675410    
## F:G:A:T     2   4.27   2.134  0.6021 0.5487344    
## F:M:A:T     2  56.53  28.264  7.9758 0.0004774 ***
## G:M:A:T     2   2.23   1.114  0.3144 0.7306023    
## F:G:M:A:T   2  14.74   7.371  2.0800 0.1278799    
## Residuals 183 648.50   3.544                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(oxygen_anova)
## 
## Call:
## lm(formula = Oxygen ~ F * G * M * A * T, data = abiotic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2844 -0.7169  0.0606  0.7525  5.7695 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.3882     0.8419  23.030  < 2e-16 ***
## F             1.0142     1.1906   0.852  0.39541    
## G             0.2158     1.2628   0.171  0.86450    
## M            -0.1926     1.1906  -0.162  0.87167    
## A            -2.6474     1.1906  -2.224  0.02740 *  
## T24          -2.1344     1.1906  -1.793  0.07467 .  
## T28          -1.1902     1.1906  -1.000  0.31879    
## F:G          -2.5408     1.7356  -1.464  0.14492    
## F:M          -1.6036     1.6837  -0.952  0.34215    
## G:M          -1.9018     1.7356  -1.096  0.27461    
## F:A           3.2560     1.7356   1.876  0.06224 .  
## G:A          -1.1100     1.7356  -0.640  0.52326    
## M:A           1.0388     1.6837   0.617  0.53803    
## F:T24         4.5275     1.7356   2.609  0.00984 ** 
## F:T28        -1.0260     1.6837  -0.609  0.54304    
## G:T24         1.9302     1.7356   1.112  0.26753    
## G:T28         0.0004     1.7356   0.000  0.99982    
## M:T24         1.4588     1.6837   0.866  0.38740    
## M:T28        -0.2414     1.6837  -0.143  0.88615    
## A:T24        -0.2418     1.6837  -0.144  0.88597    
## A:T28         1.8400     1.6837   1.093  0.27591    
## F:G:M         5.1164     2.4181   2.116  0.03571 *  
## F:G:A         3.1290     2.4544   1.275  0.20399    
## F:M:A        -2.6114     2.4181  -1.080  0.28159    
## G:M:A         3.6806     2.4181   1.522  0.12971    
## F:G:T24      -5.3471     2.4544  -2.179  0.03064 *  
## F:G:T28       0.6448     2.4181   0.267  0.79003    
## F:M:T24      -4.8777     2.4181  -2.017  0.04514 *  
## F:M:T28       0.7728     2.3812   0.325  0.74589    
## G:M:T24      -2.0615     2.5139  -0.820  0.41325    
## G:M:T28       2.1986     2.4181   0.909  0.36442    
## F:A:T24      -5.9901     2.4903  -2.405  0.01715 *  
## F:A:T28      -3.6768     2.4181  -1.521  0.13010    
## G:A:T24      -1.4251     2.4544  -0.581  0.56221    
## G:A:T28       0.9438     2.4181   0.390  0.69676    
## M:A:T24      -2.8442     2.3812  -1.194  0.23384    
## M:A:T28      -2.6162     2.3812  -1.099  0.27334    
## F:G:M:A      -6.1012     3.4197  -1.784  0.07606 .  
## F:G:M:T24     2.8517     3.5134   0.812  0.41804    
## F:G:M:T28    -4.6620     3.3937  -1.374  0.17121    
## F:G:A:T24     3.3411     3.4965   0.956  0.34056    
## F:G:A:T28    -1.5430     3.4197  -0.451  0.65237    
## F:M:A:T24     9.2637     3.4455   2.689  0.00784 ** 
## F:M:A:T28     4.5152     3.3937   1.330  0.18502    
## G:M:A:T24     0.7278     3.4881   0.209  0.83494    
## G:M:A:T28    -5.4942     3.3937  -1.619  0.10718    
## F:G:M:A:T24  -1.7131     4.9389  -0.347  0.72909    
## F:G:M:A:T28   7.5692     4.7994   1.577  0.11650    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.882 on 183 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.5535, Adjusted R-squared:  0.4388 
## F-statistic: 4.826 on 47 and 183 DF,  p-value: 6.908e-15
tc_anova <- lm(TC ~F*G*M*A*T, data = abiotic)
autoplot(oxygen_anova)

hist(residuals(oxygen_anova))

anova(oxygen_anova)
## Analysis of Variance Table
## 
## Response: Oxygen
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## F           1  24.59  24.594  6.9403 0.0091497 ** 
## G           1  21.31  21.315  6.0148 0.0151244 *  
## M           1  47.88  47.878 13.5108 0.0003115 ***
## A           1 161.67 161.674 45.6227 1.845e-10 ***
## T           2 183.84  91.919 25.9385 1.209e-10 ***
## F:G         1  12.87  12.866  3.6306 0.0582953 .  
## F:M         1   6.04   6.037  1.7036 0.1934602    
## G:M         1   1.37   1.366  0.3856 0.5353968    
## F:A         1  52.70  52.697 14.8705 0.0001593 ***
## G:A         1   4.33   4.334  1.2231 0.2702083    
## M:A         1   0.23   0.231  0.0651 0.7989686    
## F:T         2  34.02  17.012  4.8006 0.0092891 ** 
## G:T         2   9.56   4.779  1.3485 0.2622094    
## M:T         2   3.51   1.756  0.4955 0.6100538    
## A:T         2  50.91  25.454  7.1828 0.0009930 ***
## F:G:M       1  20.32  20.323  5.7349 0.0176401 *  
## F:G:A       1   8.83   8.828  2.4913 0.1162046    
## F:M:A       1   0.10   0.099  0.0278 0.8676644    
## G:M:A       1   0.00   0.001  0.0003 0.9870868    
## F:G:T       2  20.75  10.377  2.9283 0.0559984 .  
## F:M:T       2  17.88   8.939  2.5225 0.0830454 .  
## G:M:T       2   2.61   1.304  0.3681 0.6925778    
## F:A:T       2   0.31   0.155  0.0437 0.9572609    
## G:A:T       2   1.60   0.802  0.2264 0.7976449    
## M:A:T       2  18.25   9.123  2.5745 0.0789470 .  
## F:G:M:A     1  15.10  15.105  4.2624 0.0403759 *  
## F:G:M:T     2   5.41   2.705  0.7634 0.4675410    
## F:G:A:T     2   4.27   2.134  0.6021 0.5487344    
## F:M:A:T     2  56.53  28.264  7.9758 0.0004774 ***
## G:M:A:T     2   2.23   1.114  0.3144 0.7306023    
## F:G:M:A:T   2  14.74   7.371  2.0800 0.1278799    
## Residuals 183 648.50   3.544                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(oxygen_anova)
## 
## Call:
## lm(formula = Oxygen ~ F * G * M * A * T, data = abiotic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2844 -0.7169  0.0606  0.7525  5.7695 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.3882     0.8419  23.030  < 2e-16 ***
## F             1.0142     1.1906   0.852  0.39541    
## G             0.2158     1.2628   0.171  0.86450    
## M            -0.1926     1.1906  -0.162  0.87167    
## A            -2.6474     1.1906  -2.224  0.02740 *  
## T24          -2.1344     1.1906  -1.793  0.07467 .  
## T28          -1.1902     1.1906  -1.000  0.31879    
## F:G          -2.5408     1.7356  -1.464  0.14492    
## F:M          -1.6036     1.6837  -0.952  0.34215    
## G:M          -1.9018     1.7356  -1.096  0.27461    
## F:A           3.2560     1.7356   1.876  0.06224 .  
## G:A          -1.1100     1.7356  -0.640  0.52326    
## M:A           1.0388     1.6837   0.617  0.53803    
## F:T24         4.5275     1.7356   2.609  0.00984 ** 
## F:T28        -1.0260     1.6837  -0.609  0.54304    
## G:T24         1.9302     1.7356   1.112  0.26753    
## G:T28         0.0004     1.7356   0.000  0.99982    
## M:T24         1.4588     1.6837   0.866  0.38740    
## M:T28        -0.2414     1.6837  -0.143  0.88615    
## A:T24        -0.2418     1.6837  -0.144  0.88597    
## A:T28         1.8400     1.6837   1.093  0.27591    
## F:G:M         5.1164     2.4181   2.116  0.03571 *  
## F:G:A         3.1290     2.4544   1.275  0.20399    
## F:M:A        -2.6114     2.4181  -1.080  0.28159    
## G:M:A         3.6806     2.4181   1.522  0.12971    
## F:G:T24      -5.3471     2.4544  -2.179  0.03064 *  
## F:G:T28       0.6448     2.4181   0.267  0.79003    
## F:M:T24      -4.8777     2.4181  -2.017  0.04514 *  
## F:M:T28       0.7728     2.3812   0.325  0.74589    
## G:M:T24      -2.0615     2.5139  -0.820  0.41325    
## G:M:T28       2.1986     2.4181   0.909  0.36442    
## F:A:T24      -5.9901     2.4903  -2.405  0.01715 *  
## F:A:T28      -3.6768     2.4181  -1.521  0.13010    
## G:A:T24      -1.4251     2.4544  -0.581  0.56221    
## G:A:T28       0.9438     2.4181   0.390  0.69676    
## M:A:T24      -2.8442     2.3812  -1.194  0.23384    
## M:A:T28      -2.6162     2.3812  -1.099  0.27334    
## F:G:M:A      -6.1012     3.4197  -1.784  0.07606 .  
## F:G:M:T24     2.8517     3.5134   0.812  0.41804    
## F:G:M:T28    -4.6620     3.3937  -1.374  0.17121    
## F:G:A:T24     3.3411     3.4965   0.956  0.34056    
## F:G:A:T28    -1.5430     3.4197  -0.451  0.65237    
## F:M:A:T24     9.2637     3.4455   2.689  0.00784 ** 
## F:M:A:T28     4.5152     3.3937   1.330  0.18502    
## G:M:A:T24     0.7278     3.4881   0.209  0.83494    
## G:M:A:T28    -5.4942     3.3937  -1.619  0.10718    
## F:G:M:A:T24  -1.7131     4.9389  -0.347  0.72909    
## F:G:M:A:T28   7.5692     4.7994   1.577  0.11650    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.882 on 183 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.5535, Adjusted R-squared:  0.4388 
## F-statistic: 4.826 on 47 and 183 DF,  p-value: 6.908e-15
tn_anova <- lm(TN ~F*G*M*A*T, data = abiotic)
autoplot(tn_anova)

hist(residuals(tn_anova))

anova(tn_anova)
## Analysis of Variance Table
## 
## Response: TN
##            Df Sum Sq Mean Sq  F value Pr(>F)    
## F           1 459470  459470 724.3140 <2e-16 ***
## G           1     12      12   0.0195 0.8891    
## M           1     77      77   0.1212 0.7282    
## A           1      1       1   0.0021 0.9634    
## T           2   1964     982   1.5478 0.2155    
## F:G         1     12      12   0.0194 0.8893    
## F:M         1    132     132   0.2080 0.6489    
## G:M         1    312     312   0.4919 0.4840    
## F:A         1    173     173   0.2722 0.6025    
## G:A         1    285     285   0.4493 0.5035    
## M:A         1    143     143   0.2248 0.6360    
## F:T         2    992     496   0.7818 0.4591    
## G:T         2    129      65   0.1019 0.9032    
## M:T         2    250     125   0.1970 0.8214    
## A:T         2   2701    1351   2.1291 0.1219    
## F:G:M       1     59      59   0.0932 0.7605    
## F:G:A       1     54      54   0.0851 0.7708    
## F:M:A       1      2       2   0.0038 0.9509    
## G:M:A       1     56      56   0.0881 0.7669    
## F:G:T       2    573     286   0.4513 0.6375    
## F:M:T       2    139      70   0.1097 0.8962    
## G:M:T       2     17       8   0.0132 0.9869    
## F:A:T       2   1570     785   1.2376 0.2925    
## G:A:T       2    223     111   0.1756 0.8391    
## M:A:T       2     79      39   0.0622 0.9397    
## F:G:M:A     1      8       8   0.0121 0.9124    
## F:G:M:T     2    121      61   0.0957 0.9088    
## F:G:A:T     2    176      88   0.1390 0.8704    
## F:M:A:T     2    523     262   0.4126 0.6625    
## G:M:A:T     2    195      98   0.1537 0.8576    
## F:G:M:A:T   2    439     219   0.3457 0.7082    
## Residuals 183 116086     634                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(tn_anova)
## 
## Call:
## lm(formula = TN ~ F * G * M * A * T, data = abiotic)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -93.982  -0.366   0.085   1.879  61.348 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.21800   11.26368   0.019    0.985    
## F            93.77400   15.92925   5.887 1.84e-08 ***
## G             0.31700   16.89552   0.019    0.985    
## M             1.79800   15.92925   0.113    0.910    
## A             0.08000   15.92925   0.005    0.996    
## T24           0.11200   15.92925   0.007    0.994    
## T28           0.20400   15.92925   0.013    0.990    
## F:G          -6.49300   23.22067  -0.280    0.780    
## F:M          -3.72000   22.52736  -0.165    0.869    
## G:M           0.25100   23.22067   0.011    0.991    
## F:A           0.16050   23.22067   0.007    0.994    
## G:A          -0.12100   23.22067  -0.005    0.996    
## M:A          -0.11800   22.52736  -0.005    0.996    
## F:T24        -1.83400   23.22067  -0.079    0.937    
## F:T28        -2.59600   22.52736  -0.115    0.908    
## G:T24        -0.17100   23.22067  -0.007    0.994    
## G:T28         0.24500   23.22067   0.011    0.992    
## M:T24         0.10200   22.52736   0.005    0.996    
## M:T28        -0.68800   22.52736  -0.031    0.976    
## A:T24        16.93200   22.52736   0.752    0.453    
## A:T28        -0.34800   22.52736  -0.015    0.988    
## F:G:M         9.42700   32.35245   0.291    0.771    
## F:G:A         7.01250   32.83898   0.214    0.831    
## F:M:A         1.90950   32.35245   0.059    0.953    
## G:M:A         0.29300   32.35245   0.009    0.993    
## F:G:T24      10.56300   32.83898   0.322    0.748    
## F:G:T28       8.94700   32.35245   0.277    0.782    
## F:M:T24      -9.63000   32.35245  -0.298    0.766    
## F:M:T28       5.18200   31.85849   0.163    0.871    
## G:M:T24      -0.06033   33.63423  -0.002    0.999    
## G:M:T28      -0.62700   32.35245  -0.019    0.985    
## F:A:T24      -9.48250   33.31841  -0.285    0.776    
## F:A:T28       2.31150   32.35245   0.071    0.943    
## G:A:T24     -16.23950   32.83898  -0.495    0.622    
## G:A:T28      -0.23900   32.35245  -0.007    0.994    
## M:A:T24     -16.33400   31.85849  -0.513    0.609    
## M:A:T28       0.78000   31.85849   0.024    0.980    
## F:G:M:A      -9.15450   45.75328  -0.200    0.842    
## F:G:M:T24    -0.20617   47.00702  -0.004    0.997    
## F:G:M:T28   -12.27500   45.40534  -0.270    0.787    
## F:G:A:T24     8.34800   46.78157   0.178    0.859    
## F:G:A:T28   -33.16650   45.75328  -0.725    0.469    
## F:M:A:T24    26.51850   46.09859   0.575    0.566    
## F:M:A:T28   -29.10350   45.40534  -0.641    0.522    
## G:M:A:T24    16.12083   46.66843   0.345    0.730    
## G:M:A:T28     0.00900   45.40534   0.000    1.000    
## F:G:M:A:T24 -17.83783   66.07918  -0.270    0.788    
## F:G:M:A:T28  35.37850   64.21284   0.551    0.582    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.19 on 183 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.8022, Adjusted R-squared:  0.7514 
## F-statistic: 15.79 on 47 and 183 DF,  p-value: < 2.2e-16
#abiotic plots
abiotic <- abiotic %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))


stderr <- function(x, na.rm=FALSE) {
  if (na.rm) x <- na.omit(x)
  sqrt(var(x)/length(x))
}

#calculating mean of oxygen and pH 
abiotic_means <- abiotic %>%
  filter(T %in% c("20", "24", "28")) %>%
    dplyr::group_by(Treatment, Combinations, T, Stressor,Stressor_numeric) %>%
    dplyr::summarize(Oxygen_mean = mean(Oxygen),
                     pH_mean = mean(pH),
                     Oxygen_SE = stderr(Oxygen),
                     pH_SE = stderr(pH),
                     TN_mean = mean(TN),
                     TC_mean = mean(TC),
                     TN_SE = stderr(TN),
                     TC_SE = stderr(TC))
## `summarise()` has grouped output by 'Treatment', 'Combinations', 'T', 'Stressor'. You can override using the `.groups` argument.
plot_oxygen_means <- abiotic_means %>%
ggplot() +
  geom_point(aes(x=T, y=Oxygen_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=Oxygen_mean,col=Combinations,ymin=Oxygen_mean-Oxygen_SE, ymax=Oxygen_mean+Oxygen_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=T, y=Oxygen_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature") +
  theme(legend.position = "none") +
  ylab("Oxygen[%] (mean)")

plot_oxygen_numeric_means <- abiotic_means %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=Oxygen_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=Oxygen_mean,col=T,ymin=Oxygen_mean-Oxygen_SE, ymax=Oxygen_mean+Oxygen_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=Stressor_numeric, y=Oxygen_mean, col=T),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
  ylab("Oxygen[%] (mean)") +
  xlab("Number of stressors")+
  theme(legend.position = "none")

#plot oxygen
plot_oxygen_means + plot_oxygen_numeric_means

#pH
plot_pH_means <- abiotic_means %>%
ggplot() +
  geom_point(aes(x=T, y=pH_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=pH_mean,col=Combinations,ymin=pH_mean-pH_SE, ymax=pH_mean+pH_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=T, y=pH_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature") +
  ylab("pH (mean)")+
 theme(legend.position = "none")

plot_pH_means_numeric <- abiotic_means %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=pH_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=pH_mean,col=T,ymin=pH_mean-pH_SE, ymax=pH_mean+pH_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=Stressor_numeric, y=pH_mean, col=T),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
  ylab("pH (mean)")+ 
  xlab("Number of stressors")+
  theme(legend.position = "none")
  
#plot pH
plot_pH_means + plot_pH_means_numeric 

#total nitrogen

plot_TN_means <- abiotic_means %>%
ggplot() +
  geom_point(aes(x=T, y=TN_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=TN_mean,col=Combinations,ymin=TN_mean-TN_SE, ymax=TN_mean+TN_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=T, y=TN_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature") +
  ylab("TN (mean)")+
 theme(legend.position = "none")

plot_TN_means_numeric <- abiotic_means %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=TN_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=TN_mean,col=T,ymin=TN_mean-TN_SE, ymax=TN_mean+TN_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=Stressor_numeric, y=TN_mean, col=T),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
  ylab("TN (mean)")+ 
  xlab("Number of stressors")+
  theme(legend.position = "none")

#plot nitrogen

plot_TN_means +theme(legend.position = "bottom",
          legend.box = "vertical")  + plot_TN_means_numeric + theme(legend.position = "bottom",
          legend.box = "vertical")

#total carbon 

plot_TC_means <- abiotic_means %>%
ggplot() +
  geom_point(aes(x=T, y=TC_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=TC_mean,col=Combinations,ymin=TC_mean-TC_SE, ymax=TC_mean+TC_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=T, y=TC_mean, col=Combinations),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
 scale_colour_viridis_d(direction = -1)+
  xlab("Temperature") +
  ylab("TC (mean)")+
 theme(legend.position = "none")

plot_TC_means_numeric <- abiotic_means %>%
ggplot() +
  geom_point(aes(x=Stressor_numeric, y=TC_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=TC_mean,col=T,ymin=TC_mean-TC_SE, ymax=TC_mean+TC_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
  geom_line(aes(group=Combinations, x=Stressor_numeric, y=TC_mean, col=T),position = position_dodge(width = 0.4)) +
  theme_bw() + theme(panel.grid = element_blank()) +
    theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
  ylab("TC (mean)")+ 
  xlab("Number of stressors")+
  theme(legend.position = "none")
  
#plot total carbon

plot_TC_means + plot_TC_means_numeric 

final plot for figure

(plot_diversity_means +  plot_diversity_means_numeric)/
(nmds1_plot_temperature +  plot_NMDS1_means_numeric)/
(nmds1_plot_temperature_subset +  plot_NMDS1_means_numeric_subset)/
(plot_oxygen_means + plot_oxygen_numeric_means)/
(plot_pH_means +theme(legend.position = "bottom",
          legend.box = "vertical",legend.text = element_text(size=12),legend.title = element_text(size=15)) +labs(col = "Stressor combinations",shape= "") + plot_pH_means_numeric + theme(legend.position = "bottom",
          legend.box = "vertical",legend.text = element_text(size=12),legend.title = element_text(size=15))+ labs(col = "Temperature")+ labs(shape= "Temperature")) +plot_annotation(tag_levels = "a") 

#supplementary plot of all variables

(NMDS2_plot_temperature +  plot_NMDS2_means_numeric) /
(NMDS2_plot_temperature_subset +  plot_NMDS2_means_numeric_subset)/
  (plot_TN_means + plot_TN_means_numeric) /
  (plot_TC_means  +theme(legend.position = "bottom",
          legend.box = "vertical",legend.text = element_text(size=12),legend.title = element_text(size=15)) +labs(col = "Stressor combinations",shape= "") + plot_TC_means_numeric+ theme(legend.position = "bottom",
          legend.box = "vertical",legend.text = element_text(size=12),legend.title = element_text(size=15))+ labs(col = "Temperature") + labs(shape= "Temperature")) + plot_annotation(tag_levels = "a")

extract the f and t test values, make one data frame for all factors

#NMDS1
ano_NMDS1_t <- tidy(anova_NMDS1)
names(ano_NMDS1_t)[names(ano_NMDS1_t) == "statistic"] <- "t_statistic"
names(ano_NMDS1_t)[names(ano_NMDS1_t) == "p.value"] <- "t_p.value"
names(ano_NMDS1_t)[names(ano_NMDS1_t) == "term"] <- "treatment_combination"

ano_NMDS1_t <- ano_NMDS1_t %>%
  mutate(term = removeNumbers(treatment_combination))

ano_NMDS1_f <- tidy(anova(anova_NMDS1))
names(ano_NMDS1_f)[names(ano_NMDS1_f) == "statistic"] <- "f_statistic"
names(ano_NMDS1_f)[names(ano_NMDS1_f) == "p.value"] <- "f_p.value"

nmds1 <- full_join(ano_NMDS1_t, ano_NMDS1_f)
## Joining, by = "term"
nmds1 <- nmds1 %>%
  mutate(color = ifelse(nmds1$f_p.value < 0.05, estimate, NA))

nmds1$variable="NMDS1"

nmds1 <- nmds1 %>%
  filter(term != "(Intercept)", term!="Residuals")

nmds1 <- nmds1 %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))

#NMDS2
ano_NMDS2_t <- tidy(anova_NMDS2)
names(ano_NMDS2_t)[names(ano_NMDS2_t) == "statistic"] <- "t_statistic"
names(ano_NMDS2_t)[names(ano_NMDS2_t) == "p.value"] <- "t_p.value"
names(ano_NMDS2_t)[names(ano_NMDS2_t) == "term"] <- "treatment_combination"

ano_NMDS2_t <- ano_NMDS2_t %>%
  mutate(term = removeNumbers(treatment_combination))

ano_NMDS2_f <- tidy(anova(anova_NMDS2))
names(ano_NMDS2_f)[names(ano_NMDS2_f) == "statistic"] <- "f_statistic"
names(ano_NMDS2_f)[names(ano_NMDS2_f) == "p.value"] <- "f_p.value"

nmds2 <- full_join(ano_NMDS2_t, ano_NMDS2_f)
## Joining, by = "term"
nmds2 <- nmds2 %>%
  mutate(color = ifelse(nmds2$f_p.value < 0.05, estimate/1.9, NA))

nmds2$variable="NMDS2"

nmds2 <- nmds2 %>%
  filter(term != "(Intercept)", term!="Residuals")

nmds2 <- nmds2 %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
 
# NMDS1_CB_SB /1.8
ano_nmds1_mixed_t <- tidy(anova_NMDS1_mixed)
names(ano_nmds1_mixed_t )[names(ano_nmds1_mixed_t ) == "statistic"] <- "t_statistic"
names(ano_nmds1_mixed_t )[names(ano_nmds1_mixed_t ) == "p.value"] <- "t_p.value"
names(ano_nmds1_mixed_t )[names(ano_nmds1_mixed_t ) == "term"] <- "treatment_combination"

ano_nmds1_mixed_t   <- ano_nmds1_mixed_t   %>%
  mutate(term = removeNumbers(treatment_combination))

ano_nmds1_mixed_f <- tidy(anova(anova_NMDS1_mixed))
names(ano_nmds1_mixed_f )[names(ano_nmds1_mixed_f ) == "statistic"] <- "f_statistic"
names(ano_nmds1_mixed_f )[names(ano_nmds1_mixed_f ) == "p.value"] <- "f_p.value"

nmds1_mixed <- full_join(ano_nmds1_mixed_t  , ano_nmds1_mixed_f)
## Joining, by = "term"
nmds1_mixed <- nmds1_mixed %>%
  mutate(color = ifelse(nmds1_mixed$f_p.value < 0.05, estimate/1.8, NA))

nmds1_mixed$variable="NMDS1_CB_SB"

nmds1_mixed <- nmds1_mixed %>%
  filter(term != "(Intercept)", term!="Residuals")

nmds1_mixed <- nmds1_mixed %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))

#NMDS2_CB_SB

ano_nmds2_mixed_t <- tidy(anova_NMDS2_mixed)
names(ano_nmds2_mixed_t )[names(ano_nmds2_mixed_t ) == "statistic"] <- "t_statistic"
names(ano_nmds2_mixed_t )[names(ano_nmds2_mixed_t ) == "p.value"] <- "t_p.value"
names(ano_nmds2_mixed_t )[names(ano_nmds2_mixed_t ) == "term"] <- "treatment_combination"

ano_nmds2_mixed_t   <- ano_nmds2_mixed_t   %>%
  mutate(term = removeNumbers(treatment_combination))

ano_nmds2_mixed_f <- tidy(anova(anova_NMDS2_mixed))
names(ano_nmds2_mixed_f )[names(ano_nmds2_mixed_f ) == "statistic"] <- "f_statistic"
names(ano_nmds2_mixed_f )[names(ano_nmds2_mixed_f ) == "p.value"] <- "f_p.value"

nmds2_mixed <- full_join(ano_nmds2_mixed_t  , ano_nmds2_mixed_f)
## Joining, by = "term"
nmds2_mixed <- nmds2_mixed %>%
  mutate(color = ifelse(nmds2_mixed$f_p.value < 0.05, estimate, NA))

nmds2_mixed$variable="NMDS2_CB_SB"

nmds2_mixed <- nmds2_mixed %>%
  filter(term != "(Intercept)", term!="Residuals")

nmds2_mixed <- nmds2_mixed %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))

# Richness Shannon
ano_richness_t <- tidy(anova_diversity)
names(ano_richness_t)[names(ano_richness_t) == "statistic"] <- "t_statistic"
names(ano_richness_t)[names(ano_richness_t) == "p.value"] <- "t_p.value"
names(ano_richness_t)[names(ano_richness_t) == "term"] <- "treatment_combination"

ano_richness_t <-ano_richness_t %>%
  mutate(term = removeNumbers(treatment_combination))

ano_richness_f <- tidy(anova(anova_diversity))
names(ano_richness_f )[names(ano_richness_f ) == "statistic"] <- "f_statistic"
names(ano_richness_f)[names(ano_richness_f ) == "p.value"] <- "f_p.value"

richness <- full_join(ano_richness_t, ano_richness_f)
## Joining, by = "term"
richness <- richness %>%
  mutate(color = ifelse(richness$f_p.value < 0.05, estimate, NA))

richness$variable="Shannon_richness"

richness <- richness %>%
  filter(term != "(Intercept)", term!="Residuals")

richness <- richness %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))

# Tychonema Cyanobacteria /2.2

ano_tychonea_t <- tidy(anova_tychonea)
names(ano_tychonea_t)[names(ano_tychonea_t) == "statistic"] <- "t_statistic"
names(ano_tychonea_t)[names(ano_tychonea_t) == "p.value"] <- "t_p.value"
names(ano_tychonea_t)[names(ano_tychonea_t) == "term"] <- "treatment_combination"

ano_tychonea_t <-ano_tychonea_t %>%
  mutate(term = removeNumbers(treatment_combination))

ano_tychonea_f <- tidy(anova(anova_tychonea))
names(ano_tychonea_f )[names(ano_tychonea_f ) == "statistic"] <- "f_statistic"
names(ano_tychonea_f)[names(ano_tychonea_f ) == "p.value"] <- "f_p.value"

tyochonea <- full_join(ano_tychonea_t, ano_tychonea_f)
## Joining, by = "term"
tyochonea <- tyochonea %>%
  mutate(color = ifelse(tyochonea$f_p.value < 0.05, estimate/2.2, NA))

tyochonea$variable="Tychonema"

tyochonea <- tyochonea %>%
  filter(term != "(Intercept)", term!="Residuals")

tyochonea <- tyochonea %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))


# Sulfuricurvum /1.27
ano_sulfuri_t <- tidy(anova_sulfuri)
names(ano_sulfuri_t)[names(ano_sulfuri_t) == "statistic"] <- "t_statistic"
names(ano_sulfuri_t)[names(ano_sulfuri_t) == "p.value"] <- "t_p.value"
names(ano_sulfuri_t)[names(ano_sulfuri_t) == "term"] <- "treatment_combination"

ano_sulfuri_t <-ano_sulfuri_t %>%
  mutate(term = removeNumbers(treatment_combination))

ano_sulfuri_f <- tidy(anova(anova_sulfuri))
names(ano_sulfuri_f )[names(ano_sulfuri_f ) == "statistic"] <- "f_statistic"
names(ano_sulfuri_f)[names(ano_sulfuri_f ) == "p.value"] <- "f_p.value"

sulfuri <- full_join(ano_sulfuri_t, ano_sulfuri_f)
## Joining, by = "term"
sulfuri <- sulfuri %>%
  mutate(color = ifelse(sulfuri$f_p.value < 0.05, estimate/1.27, NA))

sulfuri$variable="Sulfuricurvum"

sulfuri <- sulfuri %>%
  filter(term != "(Intercept)", term!="Residuals")

sulfuri <- sulfuri %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))

#oxygen /9.26

ano_oxygen_t <- tidy(oxygen_anova)
names(ano_oxygen_t )[names(ano_oxygen_t ) == "statistic"] <- "t_statistic"
names(ano_oxygen_t )[names(ano_oxygen_t ) == "p.value"] <- "t_p.value"
names(ano_oxygen_t )[names(ano_oxygen_t ) == "term"] <- "treatment_combination"

ano_oxygen_t  <- ano_oxygen_t  %>%
  mutate(term = removeNumbers(treatment_combination))

ano_oxygen_f <- tidy(anova(oxygen_anova))
names(ano_oxygen_f)[names(ano_oxygen_f) == "statistic"] <- "f_statistic"
names(ano_oxygen_f)[names(ano_oxygen_f) == "p.value"] <- "f_p.value"

oxygen <- full_join(ano_oxygen_t , ano_oxygen_f)
## Joining, by = "term"
oxygen <- oxygen %>%
  mutate(color = ifelse(oxygen$f_p.value < 0.05, estimate/9.26, NA))

oxygen$variable="oxygen"

oxygen <- oxygen %>%
  filter(term != "(Intercept)", term!="Residuals")

oxygen <- oxygen %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))

#pH /2.62

ano_pH_t <- tidy(pH_anova)
names(ano_pH_t)[names(ano_pH_t) == "statistic"] <- "t_statistic"
names(ano_pH_t)[names(ano_pH_t) == "p.value"] <- "t_p.value"
names(ano_pH_t)[names(ano_pH_t) == "term"] <- "treatment_combination"

ano_pH_t <- ano_pH_t %>%
  mutate(term = removeNumbers(treatment_combination))

ano_pH_f <- tidy(anova(pH_anova))
names(ano_pH_f )[names(ano_pH_f) == "statistic"] <- "f_statistic"
names(ano_pH_f )[names(ano_pH_f) == "p.value"] <- "f_p.value"

pH <- full_join(ano_pH_t , ano_pH_f)
## Joining, by = "term"
pH <- pH %>%
  mutate(color = ifelse(pH$f_p.value < 0.05, estimate/2.62, NA))

pH$variable="pH"

pH <- pH %>%
  filter(term != "(Intercept)", term!="Residuals")

pH <- pH %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))

#TN/93.77

ano_tn_t <- tidy(tn_anova)
names(ano_tn_t )[names(ano_tn_t ) == "statistic"] <- "t_statistic"
names(ano_tn_t )[names(ano_tn_t ) == "p.value"] <- "t_p.value"
names(ano_tn_t )[names(ano_tn_t ) == "term"] <- "treatment_combination"

ano_tn_t  <- ano_tn_t  %>%
  mutate(term = removeNumbers(treatment_combination))

ano_tn_f <- tidy(anova(tn_anova))
names(ano_tn_f)[names(ano_tn_f) == "statistic"] <- "f_statistic"
names(ano_tn_f)[names(ano_tn_f) == "p.value"] <- "f_p.value"

tn <- full_join(ano_tn_t , ano_tn_f)
## Joining, by = "term"
tn <- tn %>%
  mutate(color = ifelse(tn$f_p.value < 0.05, estimate/93.77, NA))

tn$variable="tn"

tn <- tn %>%
  filter(term != "(Intercept)", term!="Residuals")

tn <- tn %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))


#TC /65
ano_tc_t <- tidy(tc_anova)
names(ano_tc_t )[names(ano_tc_t ) == "statistic"] <- "t_statistic"
names(ano_tc_t )[names(ano_tc_t ) == "p.value"] <- "t_p.value"
names(ano_tc_t )[names(ano_tc_t ) == "term"] <- "treatment_combination"

ano_tc_t  <- ano_tc_t  %>%
  mutate(term = removeNumbers(treatment_combination))

ano_tc_f <- tidy(anova(tc_anova))
names(ano_tc_f)[names(ano_tc_f) == "statistic"] <- "f_statistic"
names(ano_tc_f)[names(ano_tc_f) == "p.value"] <- "f_p.value"

tc <- full_join(ano_tc_t , ano_tc_f)
## Joining, by = "term"
tc <- tc %>%
  mutate(color = ifelse(tc$f_p.value < 0.05, estimate/65, NA))

tc$variable="tc"

tc <- tc %>%
  filter(term != "(Intercept)", term!="Residuals")

tc <- tc %>% 
   reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))

#all abiotic

all_abiotic <- rbind(oxygen, pH, tn, tc)


all_heat <- rbind(tyochonea,sulfuri,nmds1,nmds2,nmds1_mixed, nmds2_mixed,richness,oxygen, pH, tn, tc)

all_heat <- all_heat %>%
  reorder_levels(variable,order=c("Tychonema","Sulfuricurvum","NMDS1","NMDS2","NMDS1_CB_SB","NMDS2_CB_SB","Shannon_richness","oxygen", "pH", "tn", "tc"))

heat map of all factors

all_heat %>%
  filter(term != "(Intercept)", term!="Residuals")%>%
ggplot(aes(x = variable, y = treatment_combination, fill=color)) + 
    geom_tile() +
  scale_fill_gradient2(low="navy", mid="white", high="red", 
                       midpoint=0,breaks=c(-1,0,1),
                         na.value="lightgrey",
                           limits=c(-1.2,1.2)
                    ) +
  theme(strip.placement = "inside") +
  theme_bw()+
  theme(axis.text.x = element_text(angle = 50, hjust = 1,size=10)) +
  theme(axis.text.y = element_text(angle = 0, hjust = 1,size=10)) +
  ylab("model term [main effect or interaction]")

just visualization of F:A with all factors

all_heat %>%
  filter(treatment_combination %in% c("F","A","F:A","F:A:T24","F:A:T28"), variable != "tn",variable != "tc",color != "NA")%>%
  ggplot()+
    geom_point(aes(x=variable, y=color, shape=treatment_combination,col=treatment_combination),size=4)+
  theme_bw() +
   theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14)) +
  ylab("standardised estimate")+
   theme(axis.text.x = element_text(angle = 50, hjust = 1))+
  geom_hline(yintercept=0)

numeric analysis

slope

#oxygen
abiotic<- ps_relative@sam_data %>%
  as.tibble() %>%
  select(Treatment,Temperature, Combinations, Temperature_factor, Lable, Name,F,G,M,A, Oxygen, TN,TC, pH,Stressor_numeric, Stressor)
## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 object
abiotic <- abiotic %>%
  filter(Treatment!="ori")

abiotic <- abiotic %>%
  mutate(T2=Temperature^2,T=as.factor(Temperature)) %>%
  select(-Temperature)
  

abiotic$Oxygen<- (abiotic$Oxygen-mean(abiotic$Oxygen))/sd(abiotic$Oxygen)
oxygen_anova_numeric <- lm(Oxygen ~T+T:Stressor_numeric, data=abiotic)


oxy_slope<- tidy(oxygen_anova_numeric ,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")

oxy_slope$std.error <- NULL
oxy_slope$statistic <- NULL
oxy_slope$p.value <- NULL


oxy_slope$organisation="ecosystem"
oxy_slope$variable="oxygen"


#pH 

abiotic$pH <- (abiotic$pH-mean(abiotic$pH))/sd(abiotic$pH)
pH_anova_numeric <- lm(pH ~T+T:Stressor_numeric, data=abiotic)

pH_slope <- tidy(pH_anova_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
pH_slope$std.error <- NULL
pH_slope$statistic <- NULL
pH_slope$p.value <- NULL

pH_slope$organisation="ecosystem"
pH_slope$variable="pH"

#TN

abiotic$TN <- (abiotic$TN-mean(abiotic$TN))/sd(abiotic$TN)
tn_anova_numeric <- lm(TN ~T+T:Stressor_numeric, data=abiotic)


tn_slope <- tidy(tn_anova_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
tn_slope$std.error <- NULL
tn_slope$statistic <- NULL
tn_slope$p.value <- NULL

tn_slope$organisation="ecosystem"
tn_slope$variable="total nitrogen"

#tc 

abiotic$TC <- (abiotic$TC-mean(abiotic$TC))/sd(abiotic$TC)
tc_anova_numeric <- lm(TC ~T+T:Stressor_numeric, data=abiotic)

tc_slope <- tidy(tc_anova_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
tc_slope$std.error <- NULL
tc_slope$statistic <- NULL
tc_slope$p.value <- NULL

tc_slope$organisation="ecosystem"
tc_slope$variable="total carbon"


#NMDS1

mds_whole_res$NMDS1 <- (mds_whole_res$NMDS1-mean(mds_whole_res$NMDS1))/sd(mds_whole_res$NMDS1)
anova_NMDS1_numeric <- lm(NMDS1 ~T+T:Stressor_numeric, data=mds_whole_res)


nmds1_slope <- tidy(anova_NMDS1_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
nmds1_slope$std.error <- NULL
nmds1_slope$statistic <- NULL
nmds1_slope$p.value <- NULL

nmds1_slope$organisation="community"
nmds1_slope$variable="NMDS1"

#NMDS2

mds_whole_res$NMDS2 <- (mds_whole_res$NMDS2-mean(mds_whole_res$NMDS2))/sd(mds_whole_res$NMDS2)
anova_NMDS2_numeric <- lm(NMDS2 ~T+T:Stressor_numeric, data=mds_whole_res)



nmds2_slope <- tidy(anova_NMDS2_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
nmds2_slope$std.error <- NULL
nmds2_slope$statistic <- NULL
nmds2_slope$p.value <- NULL

nmds2_slope$organisation="community"
nmds2_slope$variable="NMDS2"


#NMDS1_subset
mds_whole_res_subset$NMDS1 <- (mds_whole_res_subset$NMDS1-mean(mds_whole_res_subset$NMDS1))/sd(mds_whole_res_subset$NMDS1)
anova_NMDS1_mixed_numeric <- lm(NMDS1 ~T+T:Stressor_numeric, data=mds_whole_res_subset)

nmds1_subset_slope <- tidy(anova_NMDS1_mixed_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
nmds1_subset_slope$std.error <- NULL
nmds1_subset_slope$statistic <- NULL
nmds1_subset_slope$p.value <- NULL


nmds1_subset_slope$organisation="community"
nmds1_subset_slope$variable="NMDS1_subpopulation"

#NMDS2_subset

mds_whole_res_subset$NMDS2 <- (mds_whole_res_subset$NMDS2-mean(mds_whole_res_subset$NMDS2))/sd(mds_whole_res_subset$NMDS2)
anova_NMDS2_mixed_numeric <- lm(NMDS2 ~T+T:Stressor_numeric, data=mds_whole_res_subset)
autoplot(anova_NMDS2_mixed_numeric)

nmds2_subset_slope <- tidy(anova_NMDS2_mixed_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
nmds2_subset_slope$std.error <- NULL
nmds2_subset_slope$statistic <- NULL
nmds2_subset_slope$p.value <- NULL


nmds2_subset_slope$organisation="community"
nmds2_subset_slope$variable="NMDS2_subpopulation"


#Shannon

div_raw$Shannon <- (div_raw$Shannon-mean(div_raw$Shannon))/sd(div_raw$Shannon)
anova_diversity_numeric <- lm(Shannon ~T+T:Stressor_numeric, data=div_raw)
autoplot(anova_diversity_numeric)

shannon_slope <- tidy(anova_diversity_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
shannon_slope$std.error <- NULL
shannon_slope$statistic <- NULL
shannon_slope$p.value <- NULL

shannon_slope$organisation="community"
shannon_slope$variable="species_richness"


#Tychonema

df_cyano$Abundance <- (df_cyano$Abundance-mean(df_cyano$Abundance))/sd(df_cyano$Abundance)
anova_tychonea_numeric <- lm(Abundance ~T+T:Stressor_numeric, data=df_cyano)
autoplot(anova_tychonea_numeric)

tychonema_slope <- tidy(anova_tychonea_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
tychonema_slope$std.error <- NULL
tychonema_slope$statistic <- NULL
tychonema_slope$p.value <- NULL

tychonema_slope$organisation="species"
tychonema_slope$variable="Tychonema"


#Sulfuri

df_sulfuri$Abundance <- (df_sulfuri$Abundance-mean(df_sulfuri$Abundance))/sd(df_sulfuri$Abundance)
anova_sulfuri_numeric <- lm(Abundance ~T+T:Stressor_numeric, data=df_sulfuri)
autoplot(anova_sulfuri_numeric)

sulfuri_slope <- tidy(anova_sulfuri_numeric,conf.int = TRUE) %>%
  filter(term!="(Intercept)", term!="T24", term !="T28")
sulfuri_slope$std.error <- NULL
sulfuri_slope$statistic <- NULL
sulfuri_slope$p.value <- NULL

sulfuri_slope$organisation="species"
sulfuri_slope$variable="Sulfuricurvum"


all <- rbind(tychonema_slope,sulfuri_slope,shannon_slope,nmds1_slope,nmds2_slope,nmds1_subset_slope, nmds2_subset_slope,oxy_slope,pH_slope,tn_slope,tc_slope)


all <- all %>%
  reorder_levels(variable,order=c("Tychonema","Sulfuricurvum","NMDS1","NMDS2","NMDS1_subpopulation","NMDS2_subpopulation","species_richness","oxygen","pH","total nitrogen","total carbon"))


plot_slope <-all %>%
ggplot(aes(x = organisation, y = estimate, colour = variable, shape=term)) +
  geom_point(size = 4,position = position_dodge(width = 0.8))  +
  geom_hline(aes(yintercept=0))+
geom_errorbar(aes(x=organisation,y=estimate,col=variable,ymin=conf.low, ymax=conf.high), width=0.1,alpha=0.4,position = position_dodge(width = 0.8))+
  ylab("Slope [number of stressors]") +
    theme_bw() +
  scale_x_discrete(limits = c("species","community","ecosystem")) +
  xlab("Level of organisation") +
  theme(axis.title.x = element_text(size = 15), axis.text.x = element_text(size=12) ,axis.title.y = element_text(size = 15),axis.text.y = element_text(size=12))